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COMOC: THREE DIMENSIONAL BOUNDARY REGION VARIANT 
THEORETICAL MANUAL AND USER'S GUIDE 

By 

A. J. Baker & S. W. Zelazny 
Bell Aerospace Company 


.SUMMARY 


The Three-Dimensional Boundary Region Variant of the COMOC 
computer program system solves the three-dimensional boundary 
region equations for flow of a viscous, heat conducting, multi- 
ple species, compressible fluid including combustion. The gov- 
erning partial differential equations are solved in physical 
variables, and allow complete two-dimensional diffusion in the 
plane transverse to the predominant direction of flow. The flow 
field may be external or confined, subsonic or supersonic, lam- 
inar and/or turbulent, and may contain up to nine or more dis- 
tinct species in frozen composition or undergoing equilibrium 
chemical reaction for a hydrogen/ oxygen/ ai r system. The program 
is equally applicable to computations in two- and three-dimen- 
sional boundary layer flows wherein diffusion in only one direc- 
tion is important. 

The COMOC computer program is based upon a finite element 
solution algorithm for the elliptic partial differential opera- 
tor in the parent equation system. It employs an explicit finite 
difference integration procedure to solve the resultant systems 
of first-order, ordinary differential equations. Boundary con- 
dition constraints on the normal flux and tangential distribution 
of each dependent variable are user-specifiable on arbitrarily 
disjoint segments of the solution domain closure. The solutions 
for each dependent variable, and all computed parameters, are 
established at node points lying on a speci f i ably non-regular 
computational lattice formed by plane t r i angu 1 at i on of the solu- 
tion domain. The numerical solution establishes the complete 
three-di mens i onal distributions of the three scalar velocity 
components, enthalpy, temperature, density, viscosity, and all 
applicable species mass fractions, as well as various integral 
flow parameters. Variable Pr.andtl number and species diffusion 
coefficient distributions may be utilized. Initial distributions 
of all dependent variables may be arbitrarily specified. 

This report documents the theoretical and mechanical 
structure of the computer program, and presents detailed guidance 
on adaptation of the code to solution of a particular problem. 
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Sample solutions are discussed for several problems, especially 
with respect to solution accuracy and speed as a function of pa- 
rameters under control of the user. Construction of the input 
data decks for sample problems is discussed. A programme r 1 s 
manual has been separately published [Ref. 1]. 


INTRODUCTION AND USER GUIDELINES 


The finite element methodology for numerical solution of 
i ni ti al -boundary value problems in continuum mechanics is under- 
going an explosive rate of growth. Formerly considered to be 
constrained to solution of problems in structural analysis, or 
other linear field problems wherein an equivalent extremum prin- 
ciple exists, the theoretical support is now sufficiently gen- 
eralized to render the method directly applicable to explicitly 
nonlinear problems, including the Na vi er-Stokes equations [Ref. 
2-4]. The COMOC computer program system is being developed to 
transmit this rapid theoretical progress (often couched in in- 
tricate mathematical formalism) into a viable and versatile nu- 
merical solution capability. As such, it must be applicable to 
diverse and complex problems in computational continuum mechanics 
while requiring minimal mathematical prowess on the part of the 
user. On the way to generation of this general purpose concept, 
several Variants of COMOC have been developed for specific prob- 
lem classes including transient thermal analysis [Ref. 5] and 
the two-dimensional Navier Stokes equations [Ref. 6], This re- 
port documents the developed Three-Dimensional Boundary Region 
(3DBR) Variant of COMOC, and describes its applicability to a 
wide range of practical two- and three-dimensional flow problems. 

The 3DBR Variant of COMOC solves the three-dimensional 
boundary region equations for flow of a viscous, heat conducting, 
multiple-species, compressible fluid including combustion. The 
governing partial differential equation system, developed in rec- 
tangular Cartesian coordinates from the parabolic Navi er-Stokes 
equations, allows complete diffusion in the plane perpendicular 
to the uniformly discernible predominant flow direction. The 
flow may be external or confined, subsonic or supersonic, laminar 
and/or turbulent, and can contain up to nine or more distinct 
species in frozen composition or undergoing equilibrium chemical 
reaction for a hydrogen/oxygen/air system. The finite element 
solution procedure marches the discretized equivalent of the 
governing equation system in the direction parallel to the pre- 
dominant flow. It numerically establishes the complete three- 
dimensional distributions of the three scalar velocity components, 
enthalpy, temperature, density, viscosity, and all applicable 
species mass fractions, as well as various integral flow param- 
eters. . No restrictions or simplifying assumptions are made for 
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the Prandtl number, and individual species diffusion coefficients 
are treated as variable parameters. Initial distributions of all 
dependent variables may be arbitrarily specified. Boundary con- 
dition constraints on the normal flux and tangential distribution 
of each dependent variable are user-specifiable on arbitrarily 
disjoint segments of the solution domain closure. The solutions 
for each dependent variable, and all computed parameters, are es- 
tablished at node points lying on a specifiably non- regular com- 
putational lattice formed by plane tri angul ation of the solution 
domain. 

All Variants of the COMOC system are built upon the macro- 
structure i 1 1 ustrated i n Fig. 1 . The Main executive routine al- 
locates core, using a variable dimensioning scheme, based upon 
the total degrees of freedom of the problem. The size of the 
largest problem that can be solved is thus limited (only) by the 
core size of the computer in use. The precise mix between number 
of dependent variables (and parameters), and fineness of the dis- 
cretization, is user-spe.ci f i abl e and widely variable. The Input 
module serves its standard function for all dependent variable, 
parameter, and geometric coordinate arrays. The Discretization 
mpdule forms the finite element discretization of the solution 
domain, and' eval uates all required finite element non-standard 
matrices and s tandard-matr ix multipliers. The Initialization 
module computes the remaining initial parametric data required 
to start the solution. The Integration Module constitutes the 
primary execution sequence of problem solution. It is based upon 
an integration algorithm for the column vector of unknowns of the 
solution, for which the discretized description is i ni ti al -val ued 
Calls to auxiliary routines for parameter evaluation, e.g. vis- 
cosity, Prandtl number, source terms, combustion parameters, etc. 
as specified functions of dependent and/or independent variables 
are governed by the Integration Module. The user has consider- 
able latitude to adapt COMOC to the specifics of his particular 
problem at this point, by directly inserting easily written sub- 
routines into COMOC to compute special forms of these parameters. 
The Output module is similarly addressed from the integration 
sequence and serves its standard function via a highly automated 
array display algorithm. COMOC can execute distinct problems in 
sequence and contains an automatic restart capability to continue 
s o 1 uti ons . 

The 3DBR Variant of COMOC, as a direct consequence of the 
expansive problem class to which it may be. addressed , is a fairly 
large and complex computer program. The vigor with which the 
potential user of a computer code attacks preparation of a data 
deck decreases exponentially (at least) with the thickness of the 
instruction manual. It is the intent of this user’s guide to, 
in a minimum amount, of space, present general guidelines for the 
use of COMOC-, describe the rudiments of the differential equation 
system being solved, briefly expose the basic mathematics 
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of the finite element algorithm and its numerical embodiment, 
and discuss sample solutions with respect to accuracy, solution, 
speed, and diversity. The standard test cases that accompany 
COMOC are discussed in terms of the physics of the solution as 
well as the description of data deck preparation. A large effort 
has been made to simplify and streamline data preparation and to 
require the user to specify an absolute minimum of non-physical 
or non-engineering input. The basic program contains a multitude 
of options which could potentially lead to confusion on the part 
of the user. Most of these have been suppressed, particularly 
in the output sub-program, with default instructions or values. 

The programmer's manual [Ref. 1] describes how they may be 
returned to an operational status. 

The following general guidelines will assist the potential 
user on adapting 30BR COMOC to a given problem. 

Solution Domain Configuration 

Most three-dimensional flow fields for which, 1) a predom- 
inant flow direction persists (i.e., no recirculation component), 

2) a prescribed pressure gradient can be established, and 3) no 
imbedded shocks occur, are amenable to analysis using 3DBR COMOC. 
This includes two- and three-dimensional boundary layer flows, 
certain two- and three-dimensional flows in environmental hydro- 
dynami cs , and free-, slot-, and boundary-jet injection config- 
urations typical of combustors. Boundary conditions can be ap- 
plied to the entire solution domain closure with local normal 
orthogonal to the direction of predominant flow. An initial dis- 
tribution (i ncl udi ng zero) of all dependent variables is needed 
to start the solution. However, a downstream outflow boundary 
condition is specifically not required. 

Variables and Parameters 

The computational variables are the three scalar components 
of velocity, stagnation enthalpy, and mass fraction of all iden- 
tifiable species. Perfect gas behavior is assumed. The present 
3DBR Variant solves the mainstream and one cross-plane velocity 
component as a boundary value problem; it employs the continuity 
equation to establish the remaining cross-plane velocity component. 
The program computes static temperature and density, and all 
thermophys i ca 1 properties may be temperature and mass fraction 
dependent. Unless overridden by a user provided subroutine, vis- 
cosity is computed from Sutherland's law. The Prandtl and Schmidt 
numbers may be variable. Equilibrium combustion of arbitrary 

mixtures of hydrogen, oxygen, and air can be established including 
local heat release and formation of MO. In the absence of dil- 
uents, this capability provides equilibrium gas behavior for air 
computations including dissociation. 
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Discretization 


The nature of the flows to which this Variant is applicable 
yields the req u i remen t for two-dimensional finite element dis- 
cretizations only. Since the continuity equation is employed to 
solve for a transverse velocity component, it is advantageous to 
have node columns oriented parallel to that coordinate. This re- 
quirement for grid regularity has been built into an automatic 
discretizer for 3DBR COMOC. Considering flow in an axial corner 
for example, see Fig. 2, it might be desired to use a finer grid 
near the walls where larger dependent variable gradients would 
exist. The user need specify (only) the desired incremental 
spacing between node columns and rows. The discretizer will au- 
tomatically triangulate the domain on these node point coordinates, 
and prepare the required geometric input data. This discretizer 
is not directly applicable to non- rectangu 1 ar domains; however, 

3DBR COMOC can accept discretizations formed manually or from 
other automated sources. 

Boundary Conditions 

Constraints can be imposed on the admissible behavior of 
each dependent variable and its normal flux, i.e., gradient, on 
all surfaces bounding the solution domain, see Fig. 2. These 
surfaces may constitute actual physical boundaries of the problem 
or be strictly mathematical. As an example of the latter, em- 
ploying symmetry planes to enclose a solution domain is particu- 
larly advantageous in terms of computer execution time and user 
input effort. The attendant vanishing gradient constraint is 
the automatic default value within the finite element solution 
algorithm, and its use does not require generation of any phantom 
cells or special node handling. Fixing the normal gradient in 
terms of the dependent variable is equally straightforward, and 
is useful for thermally porous walls or a slip wall boundary 
condition. Here again, no special cells or node handling is 
required on the part of the user. 

Input Preparation 

A concerted effort has been made to render input preparation 
minimal and in terms of physically meaningful variables and ex- 
pressions. However, should the solution to a dozen or more de- 
pendent variables be sought, the input deck can become of sub- 
stantial size. The program accepts input in the English system 
of units; it outputs non-dimensional i zing constants and solution 
parameters in several systems, and provides detailed output ar- 
rays of computed non-dimensional dependent variables. The pro- 
gram executes under automatic error control and will adjust inte- 
gration step-size to maintain an accurate and stable solution. The 
only user input required for this phase is the initial and final 
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//// ^ Boundary Condition Specification 


Figure 2 . Illustrative Finite Element Discretization 


User-Written Subroutines 


COMOC provides the user with considerable latitude for 
appending subroutines to perform specific parameter computations 
Included in this category are pressure gradient, laminar and/or 
turbulent viscosity, and Prandtl and Schmidt numbers. In all 
cases, a skeletal subroutine is filled in by the user to include 
an equation or tabular data of the parametric dependence on any 
number of independent or dependent variables in any combination. 
These subroutines are always written in terms of physical vari- 
ables with dimensions consistent with the input data. Non-dimen 
sionalization and calling sequence are controlled internally, 
and the user can obtain complete arrays of these computations 
from the output package. 




Output 


The 3DBR COMOC program contains a highly adaptive output 
subprogram. The user has considerable latitude in specifying 
output arrangements , both dimensional and non-dimensional, from 
the input deck. The output routine is adapted to compute inte- 
gral flow parameters including wall shear, Stanton number, and 
mixing efficiency. Data sets are automatically scaled and ordered 
to be geometrically similar to the physical problem for all dis- 
cretizations, both regular and non-regular. 

Computational Costs 

The computer cost associated with generating a COMOC solution 
to a given problem can be approximately estimated. CPU costs are 
basically a function of the number of dependent variables in the 
solution, the amount of output requested, out-of-core operations 
associated with restart and/or plot tape preparation', and the 
thermodynamics of the solution. The use of rather course discre- 
tizations is strongly recommended for initial evaluation of any 
problem. Employing progressively finer discretizations will gen- 
erally improve solution accuracy with a more-than-proporti onal in- 
crease in computational cost. 


NOMENCLATURE 

a boundary condition coefficient 

A species; one-dimensional matrix; area 

Ar argon 

b coefficient 

B species; two-dimensional matrix 

c coef f ic i ent 

c 0 specific heat 

C species; three-dimensional matrix 

Cf skin friction 

d differential 

D determinant 

f function of known argument 
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g function of known argument 

h static enthalpy; integration step-size; stream depth 

H stagnation enthalpy; hydrogen 

i index 

I mass defect 

AAA 

i,j,k basis vectors of rectangular coordinate system 

j summation index; coefficient 

J Jacobian 

k thermal conductivity; integration stage number; constant 

K generalized diffusion coefficient; equilibrium constant 

l differential operator; number; length 

L characteristic length; differential operator 

m number 

M Mach Number; number of finite elements 

n unit normal vector; number, nodes per element; 

dimens i onal i ty 

N nitrogen; composi tion matrix 

0 oxygen 

p pressure; predicted value 

Pr Prandtl Number 

q generalized dependent variable 

Q generalized discretized dependent variable 

r positionvector 

R domain of elliptic operator; universal gas constant 

Re Reynolds Number 

S mass source term; enthalpy boundary condition parameter 
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Sc Schmidt Number 

T temDerature 

u,U velocity 

W molecular weight 

x ^ rectangular Cartesian coordinate system 

X species mole fraction 

Y species mass fraction 

a direction cosine 

g coefficient; pressure gradient parameter 

y ratio of specific heats; turbulence i nterm i ttency 
factor 

3R closure of solution domain 

6 boundary layer thickness 

A increment 

e kinematic eddy viscosity; basis vectors; integration 

Darameter 

k coefficient 

X multiplier 

£ non-dimensional length 

p viscosity 

p density 

a integral kernel 

t integral kernel; wall shear 

<f> , $ functional 

X domain of initial value operator 

oj turbulence damping factor 

Q global solution domain 


10 



column matrix 


{ } 

[ ] square matrix 
U union 

n intersection 

Y summation 


Superscripts and Subscripts 

* approximate solution; reference state 

" derivative in the x^ direction; transf ormati on 


T 


e 

i . j » k , A 

m 

n 

o 

x 

oo 

a 


matrix transpose 
unit vector 
reference state 

constrained to solution domain closure 
effective value; local reference condition 
tensor indices 

h 

pertaining to nr subdomain (finite element) 

integration stage 

initial condition 

evaluated at x^ 

global reference condition 

species identification 

elemental species identification 
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FINITE ELEMENT SOLUTION ALGORITHM FOR THE THREE-DIMENSIONAL 

BOUNDARY REGION EQUATIONS 


The system of partial differential equations governing the 
three-dimensional boundary region flow of a compressible fluid 
is obtained from the parabolic approximation to the full Navier- 
Stokes equations. The parabolic approximation, i.e., "parabolic 
Navi er-Stokes equations," describe steady, three-dimensional flows 
wherein, 1) a predominant flow direction is uniformly discernible, 
2) in this direction (only), diffusion processes are negligible 
compared to convection, and 3) no disturbances are propagated up- 
stream antiparallel to this direction. The boundary region equa- 
tion system is obtained from parabolic Navi er-Stokes with the 
single additional assumption that a known pressure distribution 
is superimposed upon the flow field. Conversely, the approxima- 
tion may be viewed as generalization of the three-dimensional 
boundary layer equations to include diffusion processes in the 
complete two-dimensional plane of crossflow. Closure of this 
equation system requires identification of constitutive behavior. 
By employing an eddy coefficient hypothesis, the time-averaged 
turbulent flow equations appear identical to the laminar flow 
equations. Hence, the finite element development assumes a gen- 
eralized transport coefficient description, distributed as lami- 
nar or turbulent at nodes of the discretization by the user. 

The Three-Dimensional Boundary Region Equations 

In three-dimensional space, spanned by a rectangular Carte- 
sian coordinate system, identify the velocity vector 

AAA 

u i = u-ji+u 2 j + u 3 k (1) 

For development of the differential equation system, assume that 
i is aligned parallel to the predominant flow direction. Iden- 
tify a two-dimensional vector differential operator as 

( ), k = j(. ), 2 + k( ), 3 (2) 

where the comma identifies the gradient operator. Employing 
Cartesian tensor notation, with summation over 2 and 3 for re- 
peated latin subscripts, the three-dimensional boundary region 
equation system for a multiple-species, compressible, reacting 
flow takes the form 


<p u i > • i 


( pu ) 


1 ' ’1 


1 2 


0 


(3) 



pu 1 Y “ 1 

p Vl,l 

pU 1 U 3 , 1 
pu 1 H, ] 




The variables appearing in Eq. ( 3 ) - ( 7 ) are non-di me n s i ona 1 i zed 
with respect to poo, U^, Cp^, Ta>, and a length constant L, and 
have their usual interpretation in fluid mechanics. The Reynolds 
(Re), Prandtl (Pr), and Schmidt (Sc) numbers are defined with 
respect to the effective diffusion coefficient, y e , in algebraic 
combination with the laminar and turbulent contri butions as, for 
example 


ul = LL_ + 

Pr Pr Pr T 


( 8 ) 


In Eq. (8), y is the laminar viscosity, e is the kinematic eddy 
viscosity, and subscript T denotes a turbulent reference param- 
eter. The stagnation enthalpy is defined in terms of species 
staticenthalpiesas 


H = l h a Y a + u k u k (9) 

a 

The static enthalpy includes the heat of formation, h“ , of the 
species in its definition as 

h a = Tc^dT + h“ (10) 

P o 
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An equation of state is required to close the system. Assuming^ 
perfect gas behavior for each species, from Dalton's law, obtain 

P = PRT I ^ (ID 

a W a 

where R is the universal gas constant and W a is the molecular 
weight of the a-th species. 

Equilibrium combustion of hydrogen/oxygen/air systems in 
three-dimensional boundary region flow is operational in 3DBR 
COMOC. The following reactions are assumed operative. 

2H + o J H 2 0 

2H t H 2 

20 i ° 2 

H + 0 ^ OH 

N 2 + 20 £ 2N0 (12) 

The equilibrium composition of the combustion by-products is 

determined by applying the Law of Mass Action [Ref. 7] to each 
reaction defined in Eq . (12). This yields definition of a set 
of equilibrium rate constants, K, which, for the simple reaction 
nA + mB -f- S.C , are expressed in terms of species mole fraction, 
X a , as 


K E (13) 

Solution of Eq. (12) with (13), and coupled with conservation of 
total and elemental mass, yields an algebraic equation system 
for determination of the equilibrium composition of the system, 
of the form. 


[N^]{X a > = {const.} 


(14) 


O 

In Eq . (14), the elements of the matrix [N a ] account for the 
particular species mole fraction distribution, {X a }, containing 
the elemental material, e.g., 0, H, and N. 
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Finite Element Solution Algorithm 


The three-dimensional boundary region equation system, 
except for global continuity, Eq. (3), is uniformly an initial- 
boundary value problem of mathematical physics. Each of the 
partial differential equations, Eq . (4)-(7), is a special case 
of the general second-order, nonlinear partial differential 
equati on 

L ( q ) = k [ K ( q ) q , .] + f(q,q x,) - g(q,x) = 0 (15) 

5 k 11 


0 where q is a generalized dependent variable identifiable with 
each computational dependent variable. In Eq . (15), f and g are 
specified functions of their arguments, x is identified with X] 
for boundary region flows, and xi are the coordinates for which 
second order derivatives exist in the lead term. The finite 
element solution algorithm is based upon the assumption that 
L(q) is uniformly parabolic within a bounded open domain Q , , i.e., 
the lead term in Eq. (15) is uniformly elliptic within its domain 
R, with closure 3R, where 

= R x Lx 0 ,x) 06) 

and x < x < “■ Table 1 lists the functions f and g, as well 
as the appropriate parameters, for Eq. (15) identified with each 
dependent variable. 

For Eq. (15) uniformly parabolic, unique 'sol uti ons for q 
are obtained pending specification of boundary constraints on 
3R and an initial condition on RU3R. For the former, the gen- 
eral form relates the function and its normal derivative every- 
where on the closure, 3R, as 

£(q) = a^ 1 ^q(x i ,x) + a ^ 2 ^ Kq (x i ,x) , k n k - a^ 3 ^ = 0 (17) 

In Eq . (17), the a^(3T.,x) are user-specified coefficients, see 
Table 2, the superscript bar notation constrains x-j to 3R, and 
n|< is the local outward-pointing unit normal vector. For an ini- 
tial distribution, assume given throughout RU3R x x Q » 

q(x iS x 0 ) E % (x i) ( 18) 
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TABLE 2 

GENERAL BOUNDARY CONDITION STATEMENT 


Boundary Conditions 

a<’> 

a< 2 > 

a< 3 > 

No Slip at Wall 

1 

0 

0 

Slip at Wall 

+ 

1 

0 

Mass I n jecti on 

0 

1 

+ ■ 

Adiabatic Wall 

0 

T 

0 

Specified Heat Flux 

0 

1 

+ 

Temperature Dependent Flux 

+ 

1 

+ 

Symmetry Condition 

0 

1 



+ User specified as non-zero to enforce desired condition level 

Formation of the finite element solution is obtained by 
establishing the algorithm for the equation system (15)-(18). 
Straightforward theoretical development is provided by using the 
Method of Weighted Residuals (MWR) formulated on a local basis. 
Since Eq. (15) is valid throughout R, it is valid within disjoint 
interior subdomains, R m , described by (x-j ,x)eR m x txo»x) called 
"finite elements," wherein UR m = R. Form an approximate solution 
for q within R m x [xo»x)» called qjfi(xj ,x) , by expansion into a 
series solution of the form 

q*(x i ,x) = (*(x,)} T {Q(x)) m (19) 

wherein the functionals <j>|<(x-j) are members of a function set com- 
plete in R m , and the unknown expansion coefficients, Qk(x)> rep- 
resent the y-dependent values of q$(x-j,x) at specific locations 
interior to R m and on the closure, 3R m , called "nodes." Equation 
(19) is a scalar, and selection of the particular tj>k is distinct- 
ly specifiable [Ref. 8] and can be problem class dependent. 


To establish the values taken by the expansion coefficients 
in Eq. (19), require that the local error in the approximate so- 
lution to both the differential equation, L(qjfi), and the boundary 
condition statement, Jl(qjfj), for 3R m r>3R, be rendered orthogonal 
to the space of the approximation functions. Employing an un- 
known algebraic multiplier, X , the resultant equation sets can 
be combined as 


{ 4 >(x i ) }L (q*)dx - a| 


m 


<4>(x. )}A(q*)da = 0 

3R 03 R 
m 


( 20 ) 
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The number of equations (20) is identical to the number of node 
points of the finite element, R m , i.e., the number of elements, 
n, in the column matrix, {Q(x)) m > Eq • (19). 

Equation (20) forms the basic operation of the finite 
element solution. Establishment of the global solution algo- 
rithm, and determination of X, is accomplished by evaluating 
Eq. (20) in each of the M finite elements of the discretized 
solution domain, and assembly of these M x n equations into a 
global matrix system using Boolean algebra. The rank of the 
global system is less than M x n by connectivity of the finite 
element domains as well as boundary condition constraints on 3R 
where a(2), Eq . (17), vanishes identically. The lead term in 
Eq. (15) can be rearranged, using the Green-Gauss Theorem, to 
yield 


{4>(x i )}<[Kq*, k ] } dx k i){<J)(x . ) }Kq*, k n k da 


m 


3 R 


m 


- K 


(*(*,)>, kKq*,^ 


m 


( 21 ) 


For 3 RH3 R m nonvani shi ng , Eq . (21), the corresponding segment of 
the closed surface integral will cancel the boundary condition 
contribution, Eq. (20), by identifying Xa(2) with < of Eq. (15). 
The contri butions to the closed surface integral, Eq . (21), where 
3R m n3R = 0 can be made to vanish [Ref. 4]. Hence, combining Eq. 
(17)-(21), the globally assembled finite element solution algo- 
rithm for the representative partial differential equation system 
description becomes 


U 


- K 


<*}, k Kq*, k d T 
R. 


m 


{ } ( f * - g*)dr 
y v m J m ' 


m 


- K 



(1 ) * ( 3 ) \ . 

a' q* - a' '/da 
m m / 


{ 0 } 


3 R 03 R 
m 


( 22 ) 
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The rank of the global equation system, Eq. (22), is 
identical to the total number of node points on RU 8 R for which 
the dependent variable re quires solution. Equation (22) is a 
first-order, ordinary differential system, and the matrix struc- 
ture is sparse and banded. Bandwidth is a function of both se- 
lected discretization and the order of the employed approximation 
functional,'{<|>}, Eq . (19). Solution of the ordinary differential 
equation system is obtained using a finite difference numerical 
integration procedure. 

A finite element solution algorithm for the global continuity 
equation is similarly derived. Recognizing that Eq. (3) is an 
initial value problem on pU 2 as a function of X 2 » with xy and X 3 
appearing as parameters, the approximation function need span only 
the transverse coordinate direction as 


q* = {<t»(x 2 )} T {Q(x 1 ,x 3 )} m (23) 

The matrix elements Q|< are nodal values of pujj; their functional 
dependence requires so 1 uti on of Eq . (3) along lines (x],x 3 ) equal 
a constant. The solution algorithm for Eq . (3) is directly spec- 
ified as 

($}L(pu*)da =0 (24) 


where the matrix elements of {$} need not be coincidental with 
those of {<J>}, Eq . (23), and the segments R^ correspond to lines 
of ( x i , x ^ ) equal to a constant. 


-THE THREE-DIMENSIONAL BOUNDARY REGION VARIANT OF C0M0C 


The C0M0C computer program system has been established to 
embody the finite element solution algorithm for systems of equa- 
tions, Eq. (15)-(18). The computer program evaluates Eq. (22) for 
each of the appropriate dependent variables, Table 1, including 
up to nine or more species mass fractions, marches the resultant 
ordinary differential equation system downstream, and includes a 
continuity equation solver for Eq . (24). This section presents 
the theoretical aspects of these solution techniques as embodied 
in the 3DBR Variant of C0M0C. 
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Finite Element Matrix Generation 

The fundamental finite element operation is formation of 
Eq. (22) on R m and 3R m . The 3DBR Variant universally employs 
linear approximation functions, Eq. (19), for all dependent 
variables. The intrinsic finite element shapes for one- and 
two-dimensional space spanned by simplex functionals, are the 
line and triangle, Fig. 3. Accurate determination of the element 
matrices of Eq. (22) is mandatory, and involves evaluation of 
various-order moment distributions over the domain and on the 
closure of the finite element. Natural coordinate functionals, 
adapted from the area coordinates of structural mechanics [Ref. 
8], are utilized. Simplex functionals are a linearly dependent 
set of normalized functions that are orthogonal to the respective 
closure segments of the finite element domain. For an n-dimen- 
sional space, there are n + 1 simplex natural coordinate func- 
tions. Table 3 contains the implicit definition of these func- 
tions in their respective spaces. The natural coordinate 

functions vanish at all node points of the finite element except 
one where the value is unity; hence, these functions are the 
elements of the approximation functional matrix, {<{>}, Eg- (19). 
Integration of arbitrary-order products of scalar elements of 
the {<}>}, over the domain of the finite element, are analytically 
evaluated in terms of the exponent distribution, see Table 4. 

For the present case, the equation system descriptions require 
moment generation in Euclidean space spanned by a rectangular 
Cartesian basis. All computations are performed in the local 
(primed) coordinate system, Fig. 3, defined by the tensor trans- 
formation law 



One-Dimensional Space Two-Dimensional Space 


Figure 3. Intrinsic Finite Element Domains for Simplex 

Approximation Functions 
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TABLE 3 


IMPLICIT DEFINITION OF SIMPLEX NATURAL 
COORDINATE FUNCTIONS 



TABLE 4 

INTEGRALS OF NATURAL COORDINATE FUNCTION 
PRODUCTS OVER FINITE ELEMENT DOMAINS 



* D = Determinant of coefficient matrix defining the 
natural coordinate system, see Table 3. 


n = Dimensionality of the finite element space 

where ri is the position vector to the origin of the primed 
coordinate system, and the ctij are the direction cosines of the 
coordinate transformation. Tne integration kernels for two- 
dimensional space, Eq. (22), are 

dt = dx'dx' (26) 

do = dx' (27) 
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The first term in Eq . (22) is standard for all dependent 
variables. Assuming the generalized diffusion coefficient is 
distributed over the m^' 1 element as a dependent variable, obtain 


{ <D } , k Kq*, k dx 

R 


{4)}, k {K} t J{(j,}{^}I k {Q}dT 

R m 


= k{K}^{B 10}[B21 1 S ] { Q } m 


( 28 ) 


In Eq. (28) and the following, matrices with B prefixes are 
standard two-dimensional forms defined in Table 5. For Eq . (22) 
identified with each dependent variable, f* and g* universally 
contain the nonlinear convection term and the i ni xi al - val ue 
operator as dominant terms. The finite element equivalent for 
convection is 


U}pu*q* k dT = |{<()}{(t)} T {pu k } m {(j)}I k {Q} m dT 


= [B200S]{pir} m {Bll} T {Q} m (29) 

where the elements of the vector, {puF}, are nodal values of the 
planar mass flux transformed to the local coordinate systems via 

u k - a kj u j (30) 


The i ni ti al - val ue operator, which comprises the mainstream 
convection term, similarly becomes 


j^4>>pu^q* 1 dt 


(<i)H ! j)} T {pUl} m {L} T {Q} n ;dT 


= {pUl}J[B3000S]{Q}- (31) 


where the matrix elements of [B3000S] are column matrices, see 
Table 5. The superscript prime exterior to a matrix denotes an 
ordinary derivative. 
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TABLE 5 

STANDARD FINITE ELEMENT MATRIX FORMS FOR SIMPLEX 
FUNCTIONALS IN ONE- AND TWO-DIMENSIONAL SPACE 


Matrix^ 1 * 
Name 

Matrix 
Functl on 

{BIO} 

J{$)dx 

R 



[B211SJ 

{i). k {$)I k 

[B200S] 

J{+}{<f.) T do 
' R m 

[B3000S] 

]{*H*m} T dT 

R 


m 

{BID 

{*>. k 

[A200S] 

J{<»{*} T do 

3R_ 

m 

{ A 1 0 } 

J{$)do 
3 R_ 






z 1 1 

l 

li 


Matrix Eval uati on* 2 3 4 * * (3 * ’ (4) 


/ X3P3 ,V ( X3P3 -l\ 

l xTF? " ' )• j 


(1) Matrix names are a 6 digit code covering dimensionality, nonlinearity, degree of 
differentiation and special matrix properties, as [a, b, c, d, e, f] where: 

a - A, B, C for spaces of one-, two-, and three-dimensions, 
b = number of coordinate functions appearing In Integral or matrix, 
c, d, e = (0,1) Boolean counters indicating (no, yes) di fferentf atlon of each function 
e or f = S, A, A for matrix symmetric, antisymmetric or general. 

(2) Symmetric matrices are written in upper triangular form. 

(3) A ro = 1/2 (X2P2) ( X3P3) , the plane area of the triangular finite element. 

X2P2 = the x 2 prime coordinate of node 2, 

X3P3 = the x^ prime coordinate of node 3. 

(4) t ra = length of side for boundary condition (=X2P2). 








Both momentum Eq. (5) and (6) contain contributions to 
stemminq from a specified pressure distribution. For the main- 
stream momentum equation, a specified longitudinal pressure 
gradient, p , ^ , is assumed kpown; hence, 

(4, }p , dx = {B10}p»i ( x i ) (32) 


For a lateral pressure gradient, obtain 


{ 4> } P , 3 dx 


{4>>{4>}To^PZ}dx 

m 


(B10}{Bll} T {PZ> m (33) 


where the matrix elements of'{PZ) are obtained from the tensor 
transformation law, Eq . (25), as 


PZ k a k 3 P ’ 3 


(34) 


Each species continuity equation, Eq . (4), may have a source 
term. Assuming the distribution to lie over the nodes of the 
discretization, obtain 


» . 

{<j> } S a dT 


{4>H<j>} T {S a } dx = [B200S] {S a } 

III 111 


(35) 


For non-constant Prandtl and Schmidt Numbers, the energy equation, 
Eq . (7), has two source terms. An integration using a Green- 
Gauss Theorem is appropriate for both; the generated surface in- 
tegrals vanish by pairs on interior SR™ and are identically zero 
on 3 R m n 8 R for non-slip, non-porous walls. For the first term. 
Table 1 , obtai n 


U)| 


2Re 




Pr 


dx 


R 


m 


M 


2 

00 

Re 


{ ^ } ’k( I fr : / pe * u j u j, k dx 

R m 


M‘ 


- {XMUJ^PRJJBSOOOS] ^{U j } m [B2nS]{U j } n) (36) 
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In Eq. (36), the repeated subscript j is summed over all scalar 
components of the velocity vector uj . The matrix elements of 
(XMUK and { P R } m are respectively tne mth element nodal values 
of effective viscosity and the Prandtl Number function. The 
same operations repeated for the second contribution to dissi- 
pation, Eq. (7), yield 


u> 

1 / Sc-Pr\ 

Re \ScPr / 

y Z h Y . k 

t-\ 

\ / 

a 


R 


m 


dx 
5 k 


1 

Re 


UK 

R 


(Sc-Pr 
k\ScPr , 


p e *Ih“*V?*dT 


a 


m 


= - ^ {XMU}^{SC}^[B3000S][{HS a } m [B211S]{Y a } m (37) 

a ,! 

In Eq. (37), the matrix elements of (SC > m and (HS a } m are m^* 1 
element nodal values of the Schmidt Number function and the 
species static enthalpy, Eq. (10), respectively. 


The boundary condition constraint matrices, Eq . (22), are 
evaluated directly, since they are always applied on the line, 
x| equal to a constant. Using prefix A to signify a one-dimen- 
sional element operation, obtain 




UJaJqJdo 


3 R03 R 
m 


<a^ 1) [A200S]{Q} ni 


(38) 


{ } a m da 


tca^^AlO) 


(39) 


3 R f i3 R' 
m 

The A matrices are also listed in Table 5. Equations (33), (35), 
and (37)-(39) are not presently coded into 3DBR C0M0C , but are 
included here for future reference. 


Ordinary Differential Equation System Integration Algorithm 

Application of the finite element algorithm to the original 
partial differential equation has produced a large-order system 
of ordinary differential equations written on the discretized 
equivalent of the dependent variable. Several explicit numerical 
integration algorithms have been developed for ordinary differ- 
ential equations that are optimum on the multiple bases of sta- 
bility, accuracy, and required computing time, [Ref. 9, 10]. 
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However, the explicit numerical solution of stable systems of 
differential equations with large Lipschitz constants creates 
serious integration step-size restrictions. The integration 
package in COMOC contains two methods which belong to a family 
of optimally stable, 3-stage, one-step integration methods [Ref. 
11]. The operational features of this integration package, 
aside from the ease of programming using an explicit procedure, 
include being one-step (and therefore self-starting), having 
internal error control features, automatic step-size determin- 
ation, derivative evaluations required at the integration-inter- 
val end points only, and optimal stability and accuracy within 
their given structure. 


The family of numerical integration methods that are one- 
step, predictor-multiple-corrector formulas, are described by 
the equations 


pj+i 


2 

p n + l 


k- 1 
'n + 1 


a] q, 


■? 


+ hb l Pn 

+ h C b 1 PnTl + b 2 0n' ] 

, L.ruk-1 — k ~ 2 , l k ”■ 1 a ~ 1 

+ h £ b l P„ + l + b 2 


'n+l 


a k „ 
a, q, 


+ h[b!| + b$ q'] 


k 

’2 -’n- 


(40) 


Two members belonging to the 3-stage family are operational in 
COMOC. Both methods are first-order accurate, i.e., their as- 
sociated truncation error is of order h^, where h is integration 
s tep- s i ze , and they represent optimally stable methods within 
the collection of first-order accurate methods. The coefficients 
in Eq. (40) for these two methods are listed in Table 6. Method 
1 enjoys a large absolute stability interval, v/h i 1 e method 2 has 
an extended relative stability interval. Both options in the 
integration package attempt to extremize integration step size 
automatically, based upon internal error control. The estimation 
of relative truncation error for both methods is of the form 


RTE 


h_ I p n + l ” p n I 

6 L ivii 


(41 ) 


where the parameter, 3, equals 3 and 6, respecti vely , for method 
1 and 2. Equation (41) is utilized within the integration pack- 
age to evaluate the relative truncation error associated with . 
using the given integration step size, h, to estimate the (n+l) s 
value of the dependent variable. If the computed error is less 
than the user-supplied acceptable limit, the (n+l) st estimate 
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TABLE 6 


COEFFICIENTS IN INTEGRATION ALGORITHM 


FOR TWO, ONE 

-STEP, THREE-STAGE 

METHODS 

Coef f ic i en t 

Method 1 

Method 2 

a! 

1 .000000 

1 .000000 

■r 

1 . 000000 

1 .000000 


1 . 000000 

1 .000000 


1 .000000 

1 . 000000 

' 

»? 

0.037037 

0.074074 

»i 

0.962963 

0.925926 

. s 


0. 148148 

I 

0.296296 | 


0.851852 



0.703704 


for the dependent variable is accepted. Dependent upon computed 
solution behavior, the integration algorithm will selectively 
seek to increase h, by some fraction, before proceeding to the 
next integration computation. In this fashion, the package con- 
sistently seeks to increase step size (hence decrease solution 
computation time). If at any point the computed relative error 
exceeds the limit, the current predicted values for the depend- 
ent variable are discarded, a smaller step size selected, and 
the operations of Eq. (40) repeated until an acceptable error is 
measured. 

Continuity Equation Solver 

Since an explicit integration algorithm is utilized for 
solution of Eq. (22), solution of Eq. (24) for transverse mass 
flux pu2 is required only after all other dependent variable 
distributions have been obtained on the plane X] = constant. 
Establishment of (pu$) ,3 is direct since the nodal distribution 
of {pU3} is known. However, an evaluation of ( p u ^ ) , 1 is re- 
quired, since no streamwise derivatives of a dependent variable 
can be formed before the distribution of all variables is known 
in a plane. In the discretized solution, the actual requirement 
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is to establish { p U 1 } " ; the following second-order accurate 
finite difference formula for the derivative at the end point 
of two panels of data of dissimilar length is employed. 


{pUl} 


n + 1 


h n tl n + l (h n + h n + l ) 


V 2h n+ 1 + h n><PUl > n+1 


- < h „ + h n+ 1 ) 2 ipUT + h 2 +){ pUl} n ., 


(42) 


In Eq. (42), h + ; and h are the x-j integration step-sizes, 
respectively, Between trie current X] station, x n + i , and the 
previous two stations. 


An analytic expression is then established for the X 2 
distributions of mass flux derivatives, with X 3 as a parameter 
and on a nodal basis, as 

• 

( P U1 ) ' = l a U (Xq)Xp 

k= 0 K J * 

(pU3) M = l b (x )x k (43) 

J k =0 K J * 

x. L_ 

using an n L order running-smoothing polynomial generator over 
appropriate sequential panels of data. Using a unit step for 
the weighting function, 0 , Eq . (24) then takes the form 


l (pu 2 ) , 2 ' Jo^k^J + b k { X 3 ^ x 2 


dx, 


(44) 


m 


Since all terms in Eq. (44) are integrals of perfect differen- 
tials, the solution for the increment in transverse mass flux 
over an interval Ax^ is directly obtained as 

n x k+1 

A(pu£) = l Ca k (x 3 ) + b k (x 3 )] - k ~ +1 (45) 

k — 0 

Repeating Eq. (45) along each node column completes determination 
of pu* at the nodes of the transverse plane. 

Computation of Equilibrium Composition and Thermodynamic 
Properties of Hydrogen/Oxygen/Air Mixtures 

The 3DBR Variant of COMOC can compute three-dimensional 
frozen flow mixing of arbitrary gas mixtures, as well as the 
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equilibrium combustion of hydrogen/oxygen/air systems. For the 
latter, the NASA computer code GAS (see Ref. 12) has been made 
operational within COMOC after cons i derabl e modi fi cat i on to ren- 
der it compatible with a marching-type solution with multiple 
nodes, hence solutions. The equilibrium composition and thermo- 
dynamic properties of hydrogen/oxygen/air mixtures are evaluated 
as a function of temperature and pressure; relative concentra- 
tions of the elements, H ^ > 02’ ^2‘ anc * are a ^ so determined. 

The species considered are H 2 O, U 2 > H 2 , N 2 , Ar, 0, H, NO, and OH. 
Since all thermophys i ca 1 properties are temperature dependent, 
stagnation eVthalpy is typically not known a priori; consequently, 
initialization is based upon a user input total temperature dis- 
tribution. As a function of input pressure at initialization and 
the built-in tables of thermodynamic data, distributions of static 
temperature, frozen specific heat, and stagnation enthalpy cor- 
responding to input total temperature are determined using an it- 
eration algorithm based upon the method of false position. All 
solutions following i ni ti al i zati on are based upon iteration to 
equilibrium composition using computed nodal static temperature 
as the convergence parameter. The iteration on temperature is 
assumed to have converged when the difference between successive 
iterates is less than 0.1 percent. 

After convergence to a static temperature, the equilibrium 
constants for chemical reaction are calculated from the Gibbs' 
function. Composition is then determined using a modified Newton- 
Raphson iterative procedure for solution of a system of nonlinear 
algebraic equations. Once the nodal species equilibrium (or fro- 
zen) composition is determined, enthalpy, entropy, molecular 
weight, and specific heat are calculated for mixtures of ideal 
gases in terms of the computed species mole fractions, X a , as 


Molecular Weight: W 

Specific Heat: c^ 

Static Enthalpy: h 

Entropy: S 

Mass Fraction: Y a 

Gas Constant: y 


a,, a 


l x a w 

a 




a ua 


s l x h 


a 


a 


a 


,X a W a /W 


iJL 


C p - R/W 


Inp - In X° 


(46) 

(47) 

(48) 

(49) 

(50) 

(51) 
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The equilibrium model is based upon the conservation 
properties of the sum of mole fractions, and the constancy of 
the atomic number density ratios of argon/nitrogen, nitrogen/ 
oxygen, and hydrogen/oxygen. The five chemical reactions con- 
sidered are 

2H + 0 ^ H 2 0 

2H J H 2 

20 t 0 2 

H + 0 t OH 

N 2 + 20 t 2N0 (52) 


Applying the Law of Mass Action [Ref. 7] to each reaction in 
Eq. (52), the following system of nonlinear algebraic equations 
relating species mole fractions, X^ a ', is obtained. 


X 

X 

X 

X 

X 


(4) = 

K-| p^ X ^ 1 ^ 

(3) = 

k 2 p 


(6) = 

k 3 p 


(2) = 

K 4 pX* 5 M 

(7) = 

K 5 P ,/ 2 [x< 


2 X (5) 


( 8 ) 


1/2 



(53) 


In Eq. (53), K. is the equilibrium constant for the i ^ reaction, 
which is ,a function of temperature only, and p is the static 
pressure. The numbering scheme for species identification is 
listed in Table 7. 


TABLE 7 

SPECIES IDENTIFICATION FOR REACTING HYDROGEN/OXYGEN/ AIR SYSTEMS 


Number 

1 

2 

3 

4 

5 

6 

7 

8 

9 1 

Chemical Species 

H 

OH 

H 2 

h 2 o 

0 

°2 

NO 

n 2 

Arj 
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Equations for the conservation of total mass, and the individual 
atomic species H, 0, N, and A r, may be expressed in terms of 
known constants by the matrix equation. 


where 


[N p ] 
L a 


and 


[N^] { X a } 
a 


1 

1 

0 

0 

0 


1 

1 

1 

0 

0 


1 

2 

0 

0 

0 


{const. } 


{const. } 


1 

2 

1 

0 

0 


1 

0 

1 

0 

0 


1 

0 

2 

0 

0 


1 

0 

1 

1 

0 


1 . 0 

c< 2 > 

c< 3 > 

.(4) 


The specific values of the constants c 
the initial composition. 


( 6 ) 


1 

0 

0 

2 

0 


(54) 


(55) 


(56) 

are determined from 


Of the several possible choices, the computed composition 
is based upon solution of the nonlinear equilibrium equations 
for mole fraction of hydrogen, atomic oxygen, and the square root 
of molecular nitrogen. The resultant nonlinear equation system 
requiring solution is 


[f i (X a )]{X a > = {0} 


(57) 


The Newton- Raphson iteration algorithm assumes, given a set of 
trial values, X“ , determination of a new set of values Xn + 1 , 
separated from the initial estimate by AX°, by differentiating 
Eq . ( 57) to yi el d . 

[J(f i ,X a )]UX“} = - {f i (X«)} (58) 


In Eq. (58), the Jacobian contains elements,, j, , , 
numerically as 



3X“ 


determi ned 


(59) 
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The n+l st estimate of X a is accepted as the solution to Eq. (58) 
when 

] l f j ( Ci>l 1 e (60) 

J = 1 J 

_ 5 

where e is a prescribed small parameter, usually 10 . A maximum 

of thirty iterations are allowed for the solution of Eq . (57)- 
(60) to converge within e. In only a few cases has non-conver- 
gence occurred, always within a few degrees of the threshold tem- 
perature for dissociation. For these initially divergent solu- 
tions, the equations are resolved assuming that dissociation is 
negligible, i.e., the mole fractions of H, 0, OH, and NO are neg- 
ligibly small in comparison to H^, , N 2 , and H 2 O. 

ILLUSTRATIVE SOLUTIONS 


.The 3DBR Variant of C0M0C has established solutions for 
several two- and three-dimensional boundary region flows covering 
a wide range of Mach and Reynolds numbers. Several are discussed 
to illustrate the various features of solution. The data decks 
for two of these cases are presented in, the next section, and 
come as standard test cases with the program. 

Constant Density Flow Fields 

Because of its basic simplicity, the two-dimensional, iso- 
energetic, laminar boundary layer flow of a fluid at small Mach 
number (M < 0.3) provides an excellent check case for evaluating 
the essential performance features of the finite element algo- 
rithm for Eq. (15)- (17). Only one dependent variable (u-|) need 
be integrated numerically, along with solution of the continuity 
equation for Up. However, C0M0C assumes all flows are three- 
dimensional and compressible with temperature-dependent thermo- 
physical properties. The two-dimensi onal i ty is readily obtained 
by specifying only one column of elements, see Fig. 4, and en- 
forcing the vanishing normal gradient (q, n s 0) boundary condi- 
tion on the lateral segments of 3R. This is particularly simple 
since vanishing gradient is the automatic default value intrinsic 
to the finite element algorithm. The discretization may be ex- 
tended beyond the boundary layer thickness, 6 X , so that vanishing 
normal gradient may be applied along the freestream segment of 
3R as well. The slope of the diagonals of the discretization, 
Fig. 4, bears little impact on solution accuracy for two-dimen- 
sional problems. Dependent upon initial conditions and/or other 
perturbations placed into the solution, the computed variable 
distributions along each node column may differ slightly. These 
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small differences typically damp out as the solution proceeds, 
and interchanging diagonal slopes simply exchanges the location 
of these differences if they exist. The element aspect ratio 
is similarly unimportant for two-dimensional problems. For iso- 
energetic flows, COMOC contains a simple but q u i c k - ru nn i ng sub- 
routine to compute density and static temperature for two-compo- 
nent perfect gas mixtures with a specified pressure distribution. 
Unless overridden by a user-supplied subroutine, viscosity is 
computed from Sutherland's law. 



Figure 4. Finite Element Discretization for Two-Dimensional 

Boundary Layer Flow 

Assessment of computed solution accuracy and convergence 
with discretization can be obtained for the vanishing pressure 
gradient case by comparison to the B1 as i us s i mi 1 ari ty solution 
[Ref. .13]. The test case corresponds to air at atmospheric 
pressure and M = 0.272. The boundary layer thickness at the 
initial station is 6 0 = 0.0011 m. [= 0.11(-2)) and the unit 
Reynolds number is Re x = 0.63(7)/m. Two uniform discretizations 
were employed, the first spanning S 0 with only four (4) finite 
elements, while the second doubled that number. The numerical 
computations were initialized with the similar solution profiles 
at a station downstream from the leading edge at x-|/6o = 278. 

The solution was continued downstream to.x-|/6 0 = 2560. Shown 
in Fig. 5 is a comparison of the Bias i us solution to the com- 
puted velocity profiles, U] and U 2 > obtained for the coarser 
discretization. For reference purposes, the initial profiles 
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are also shown with the finite element node locations superim- 
posed. The longitudinal flow initially contained within 6 0 has 
been retarded by a factor of 2 to 3 throughout, and agreement 
between the computed and Blasius velocity profiles is excellent. 
The computed skin friction and displacement thickness distribu- 
tions are shown in Fig. 6. The computations using the coarser 
discretization slightly underpredict skin friction and over 
estimate displacement thickness. Doubling the discretization 
(to 8 elements lying within 5 0 ) noticeably improves computed 
solution agreement with the Blasius solutions. Fig. 6. Figure 7 
presents actual percent inaccuracy in the computed solutions for 
skin friction and displacement thickness. The influence of the 
coarse discretization is most noticeable in 6; however, the error 
rapidly decreases as the boundary layer grows into the discreti- 
zation, which corresponds essentially to grid refinement. As a 
function of discretization, computed skin friction, C^, converges 
approximately proportional to the square of refinement. This 
agrees exactly with the convergence rate predicted theoretically 
for the parent diffusion equation, neglecting convection, using 
linear finite element approximation functionals [Ref. 14]. The 
uniformily small inaccuracies in computed skin friction for both 
discretizations indicate that solutions, adequate for certain 
engineering approximations, can be obtained using finite element 
discretizations that appear rather coarse in comparison to con- 
ventional experience. 



Figure 5. Computed Velocity Distributions, M = 0.272, 

Re = 0 . 63( 7 ) /m 

A 
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The increased accuracy of the finer discretization solutions 
is obtained with a measurable increase in computer CPU time. 

Shown in Fig. 8 are the distributions of integration step-size, 

AX 1 , automatically established by COMOC as within the stability 
interval of the Method 1 integration algorithm for this boundary 
layer problem. The periodic stepping feature is illustrated for 
the four element solution; the step-size of the eight element 
solution remains essentially constant at about one-fourth the 
average of the former. Since the finer discretization contains 
twice as many elements as well, an approximate eight-fold in- 
crease in computer time would be anticipated. The actual in- 
crease was by a factor of 6.9, see Fig. 8. Computer CPU is also 
a function of the number of dependent variables in the solution, 
as well as user utilization of the I/O features of COMOC. An 
approximate formula for estimating execution time in seconds is 

CPU = Ajc. (61) 

i 

where x is the ratio of CP speed of the computer in use to that 
of the IBM 360/65 (e.g., X * 0.2 for the CDC 6600), and 

c-j 5 0.011 x (No. elements) x (No. passes) 

x (No. dependent variables + 0.165) 

c 2 = 0.8 x (No. outputs) x (No. pages per output) 

c^ = 1.25 x (No. outputs if Restart tape is written) 

In the expression for c-j , the number of passes is an output 
parameter from COMOC that is somewhat greater than three times 
the number of integration steps in the solution. The factor 
0.165 accounts for simple thermodynamic, viscosity, and other 
parameter evaluations. The contribution from Co occurs only 
when writing a restart tape, and the coefficient may vary with 
different computers. 

A second illustrative solution for 3DBR COMOC in the low 
speed category comes from environmental hydrodynamics. Analyti- 
cal prediction of thermal and/or waste water dispersion into 
waterways could help understand the important mechanisms for 
turbulent transport phenomena. Such a capability could also 
circumvent some of the detailed laboratory experimentation now 
required for certification of engineering projects. The example 
[Ref. 15] corresponds to turbulent dispersion of ejectant from a 
submerged waste water outfall. The solution domain is non-reg- 
ular and corresponds to the measured cross-sectional depth dis- 
tribution of a natural stream [Ref. 16] with a span of 48 m. 

(160 ft), see Fig. 9. The initial longitudinal velocity dis- 
tribution was established from the measured isovel distribution. 
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Figure 8. Integration Step-Size Distribution, 


M = 0.272, Re x = 0.63(7)/m 

Velocity m/s 



Figure 9. Cross-Section of a Natural Stream 

Showing Measured Isovels, [Ref. 16] 
by interpolation at the nodes of a 468 finite element discreti- 
zation of the cross-section, see Fig. 10. Since the solution 
domain is quite non-regular, the automatic discretizer in COMOC 
was not applicable. The waste water ejector was assumed located 
in the deepest section of the river as shown in Fig. 10. 



Figure 10. Finite Element Discretization of Stream Cross-Section 
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The flow field was assumed isoenergetic and without cross- 
flow. Hence, a marching type solution was required for main- 
stream velocity and a single species mass fraction. The total 
and static temperatures for this case are identical; thus, a 
large input pressure was employed to coerce the perfect gas sub- 
routine in COMOC to compute a uniform density distribution cor- 
responding to that of water. Closure of the governing equation 
system was obtained by specifying a turbulent viscosity law, and 
providing a user-written subroutine to override Sutherland's Law. 
A tensor turbulence law was assumed applicable [Ref. 17, 18]; 
in the plane of the finite element discretization, the eddy vis- 
cosity coefficients in the vertical (X2) and transverse (X3) 
coordinate directions were assumed given as 

y® 2 = k 2 U*h (62) 

y® 3 = k 3 U*h (63) 

where local depth of water is given by h, k 2 , and k 3 are empir- 
ical constants, and U* is the friction velocity defined as 
U* = /x/p. The approach of Patankar and Spalding [Ref. 19] was 
employed to evaluate wall shear, x, as. 

x = K 2 pU 2 [R~ ] - 0 . 1 56R~° ' 45 + 0.08723R' 0,3 

XX X 

+ 0 . 0371 3R” 0 ’ 1 8 ] (64) 

X 

0 

where K is an empirical constant (set equal to 0.435,), R x =, RK 
where R is a local Reynolds Number defined as R = ptJi^/y, U is 
local longitudinal velocity near its extremum, and X 2 j,s a rep- 
resentative length scale. For the present case, both U and X 2 
were obtained directly from the solution for the detailed veloc- 
ity profiles. A study was performed, see Table 8, to measure 
the sensitivity of the computed pollutant distribution to the 
constants in the eddy viscosity law, Eq. (62)-(63), as well as 
the tensor character. The base line case corresponds to use 
of mean depth averages for the coefficients, confirmed experi- 
mentally .to capture the essential parabol i c character of the 
measured distributions, Case II. Case III corresponds to a 
scalar eddy viscosity equal to the magnitude of the Case I 
tensor expression. 

The results obtained from this type of study are summarized 
in Fig. 11, which presents predicted mass fraction contours of 
the contaminant, for the three cases, at a station 9.6 m. down- 
stream of injection. Comparing the results of Cases I and II, 
the neglect of the parabolic distribution in the vertical mixing 
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TABLE 8 • 

PARAMETERS IN POLLUTANT DISPERSION STUDY 



k 2 

k 3 

Comments 

I 

0.067 

0.23 

Base line case [Ref. 16] 

1 1 

0. 36U-S 2 ) 

0.23 

Vertical parabolic distribution 
(£=nondimensional local depth) 

III 

0.24 

0.24 

Scalar of equal magnitude to 
Case I 



Figure 11. Predicted Mass Fraction Contours at 9.6 m Downstream 
of Injection, Three Diffusion Models 

coefficients is confirmed to be a reasonable assumption at this 
distance downstream. However, for these conditions, the omission 
of the tensorial character of the dispersion, coefficient. Case 
III, is quite measurable. Comparing Cases III and II, the larger 
vertical ( k 2 ) coefficient has allowed the 3% contour to break to 
the surface of the river. An overall larger diffusion has also 
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occurred, although, as expected, the lateral extent of the dis- 
tributions is considerably less affected. Although no experi- 
mental data are avai 1 abl e to directly confirm these results, these 
predictions amply illustrate the potential to examine trends and 
isolate key features, while capturing the important geometric 
non- regu 1 ari ti es and differential equation non-linearities so 
important to the physics of the problem. The numerical procedure 
is readily adaptable to relocation of the ejector and alteration 
of its geometry, see Fig. 12 for example. For all cases, inte- 
gration was continued downstream a distance of approximately 90 m; 
at this point the maximum ejectant concentrati on had decreased to 


M/W/i 


Initial 100% Contour 



Figure 12. Predicted Mass Fraction Contours at 9.6 m Downstream 

of Interface Injection 


15% ±1% dependent upon the viscosity law. A typical execution 
time on the IBM 360/65 was 875 s including about 145 s to produce 
an inch of output. The predicted value using Eq. (61) is 900 s. 

(It should be noted that, although tensor turbulent transport 
properties can be utilized in 3DBR, program modification, beyond 
the scope of the casual user, is required for the present Variant.) 

Compressible Flow Fields 

A standard check case for 3DBR COMOC corresponds to a nominal 
Mach 5. laminar, two-dimensional boundary layer flow over an adia- 
batic wall in a favorable pressure gradient. A similarity solu- 
tion [Ref. 20], as well as finite difference procedures, can be 
utilized to evaluate accuracy and consistency of solution trends 
for the detailed coupling of the mechanics and thermodynamics of 
this solution. The discretization and boundary condition specifi- 
cations are essentially identical to those of the small Mach num- 
ber boundary layer solution. For constant specific heat, which' 
corresponds to the similarity solution, the simple thermodynamic 
subroutine can be used to compute local static temperature and 
density for isoenergetic flow. Only the equation for u] need be 
integrated downstream, coupled with solution for U2 using the 
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continuity equation solver. For non-isoenergetic flow or variable 
specific heat, stagnation enthalpy (H) must be added to the inte- 
grated solution vector, and the local solutions for density, 
static temperature, and specific heat are obtained from the com- 
bustion subroutine. In this instance, the air composition must 
be initialized as well, although the oxygen and nitrogen elemental 
species mass fractions need not be integrated since the flow field 
composition is homogeneous and constant. For either thermodynamic 
procedure, the stagnation enthalpy is initialized from an input 

(constant) total temperature distribution. The initial U] pro- 
files are established from the similar solution for $ = 0.5 and 
S = 0 [Ref. 20]. Sutherland's law is employed to compute 
viscosity. 

The standard test case is initialized at xi = 0.03 m down- 
stream from the surface leading edge. The boundary layer thick- 
ness at this station is 6 q = 0.0039 m, the local Mach number is 
M = 3.77, the unit Reynolds number is Re x = , 83(5)/m, and the 
adiabatic wall temperature is T w = 1000°K (1800°R). Shown in 
Fig. 13 are the COMOC computed skin friction, freestream Mach 
number, and boundary layer thickness distributions for the con- 
stant specific heat case. These were obtained using two uniform 
finite element discretizations correspond!' ng to 4 and 8 elements 
spanning the initial boundary layer thickness. The input static 

pressure distribution, p e (x-]), is also presented for reference, 
and the boundary layer thickness has increased greather than 
four-fold within the solution domain. Only small differences, 
on the order of about 2%, exist between the two solutions, with 
the finer discretization producing a slightly larger skin fric- 
tion and smaller freestream Mach number. Superimposed in Fig. 

13, for. comparison purposes, are the results for the similar 
solution [Ref. 20], and a 20 zone finite difference solution ob- 
tained using the von Mises coordinate transformation [Ref. 19]. 
Agreement among the four solutions is excellent (within 2 %) for 
skin friction. The similar solution for M e lies between the 
COMOC and finite difference solutions, and agreement is within 
±3%. Shown in Fig. 14 are computed velocity profiles at X]/6 0 = 
22.7, which is about mid-way through the standard test solution 
domain. Shown for reference is the initial u-| profile obtained 
from the similar solution [Ref. 20], with the node locations of 
the 4 element discretization superimposed. Both COMOC solutions 
produce u] distributions that are slightly more concave upward 
in the mid-region in comparison to the similarity or finite 
difference solutions. The finer discretization COMOC solution 
lies closer to the similarity solution in the region where the 
two finite element solutions differ. The finite difference 
solution lies appreciably below both the COMOC and similarity 
solutions near freestream. The COMOC computed transverse veloc- 
ities are also shown in Fig. 1.4; only slight differences between 
the two discretization solutions, are apparent. The trends of 
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Mach Number - M e 
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Figure 13. Computed Supersonic Boundary Layer Parameters 
M = 5 , Re = . 83 ( 5 ) /m , 8 = 0.5 
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Figure 14. Computed Supersonic Boundary Layer Velocity 
M = 5, Re = . 83 ( 5 ) /m , 8 = 0.5 

A 


43 








the COMOC solutions are in excellent agreement with the estab- 
lished procedures; unfortunately, since each method of solution 
is distinctly numerical, no absolute accuracy assessment is 
established, as was possible for the constant density boundary 
layer check case. 

The computation of transverse velocity warrants additional 
comment. COMOC will accept, but does not require (since one is 
rarely available), an input distribution for U 2 at the initial 
station. For the discussed Mach 5 solutions, the initial U2 
distribution. Fig. 14, is self-determined by withholding its 
computation until the equation had been integrated forward 
a few stations. (This is mandatory, even if an initial U 2 
distribution is input, since several data stations are required 
to allow evaluation of Eq . (42) for (pu-|)'.) Computation of U 2 
is then initiated and it rapidly becomes consistent with the 
computed u] distributions. Solution is terminated after a few 
more steps downstream, and the computed nodal U 2 distribution 
trend with longitudinal distance is back extrapolated to estimate 
an initial distribution. Only one or two iterations of this ty be 
are typically required to establish a consistent U2 distribution. 
Starting with a zero initial distribution is probably the most 
convenient choice for analysis of engineering problems, wherein 
detailed initial accuracy is not of primary importance. 

Solution speed and accuracy have been illustrated to depend 
directly upon discretization. Solution speed is also a direct 
function of a user input control parameter (e) that places some- 
what flexible bounds on acceptable relative truncation error, 

Eq. (41), during downstream integration. A value of e of 1 0 ( - 4 ) 
has been found by experimentation to be generally adequate to 
maintain solution stability, hence accuracy. Using smaller values 
of e sharply increases computation time without producing the 
attendant increase in solution accuracy that a finer discretiza- 
tion would yield. However, increasing e up to 1 0 ( - 2 ) can produce 
measurable cost savings with a good probability of only a marginal 
decrease in solution stability and accuracy (upon rare occasion, 
the solution may actually go unstable and blow-up). Shown in 
Fig. 15 are computed integration step-size histories for the 
Mach 5 test case for several discretizations and values of e. 

In all cases, the computed solutions were of comparable accuracy 
and uniformly acceptable. The automatic stepping feature is 
illustrated, and the predicted necessary sharp decreases in Ax-j 
appears independent of e. For the two, 8-element discretization 
tests, increasing e to 1 0 ( - 2 ) about doubled the allowed step size 
and hence decreased CPU by almost 50%. As expected, a four-fold 
increase in CPU accrued for twice the discretization, for e con- 
stant at 1 0 ( - 4 ) , with integration step-size about half that of 
the baseline 8-element case. The coarsest 4-element discretiza- 
tion test recorded the largest integration step distribution with 
a remarkably smaller CPU. However, this solution appeared on the 
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Figure 15. Integration Step Size Distribution, 

M = 5, Re x = . 83 ( 5 )/m , 8 = 0.5 

ragged edge of instability, as evidenced primarily by sometimes 
anomalous behavior in transverse velocity, u^. From experimenta- 
tion, incipient overall solution instability is almost always 
foretold by erratic behavior of the computed U 2 distribution, 
since it immediately reflects anomalies in the streamwise deriv- 
ative distribution for ui , see Eq. (45). The user may typically 
feel confident that solutions generated using e>l 0 ( -4 ) are ac- 
ceptable provided the computed U2 behavior is smoothly consistent 
However, it is worthwhile to substantiate this assumption, by a 
shorter run at smaller e and/or finer discretization, if it can 
be afforded 
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Longitudinal Velocity - Uj/Ug (x i ) 


As a final evaluation, the nominal Mach 5 test case was 
repeated for temperature-dependent specific heat. This necessi- 
tates addition of stagnation enthalpy, H, to the dependent 
variable vector and addressing the thermodynamic package in the 
combustion subroutine. Shown in Fig. 16 is a comparison between 



Transverse Coordinate - x 2 /O o 

Figure 16. Computed Boundary Layer Velocity Profiles, 

M = 5, Re = . 83 ( 5 )/m , $ = 0.5 

A 

computed longitudinal velocity distributions at X]/6 Q = 22.7 for 
a uniform adiabatic wall temperature of 1000°K. Maximum differ- 
ences of about 4% accrue due to temperature dependent specific 
heat; the computed increase in skin friction is about 10%. 
Besides the increased solution cost stemming from employing H as 
a dependent variable, an additional expenditure accrues from 
addressing the more comprehensive thermodynamics package for 
variable thermophysical properties. Referring to Eq. (61), a 
fourth contribution to CPU for non-combustion utilization of 
the combustion subroutine is 


C„ = 0.004 x (No. Passes) x (No. Nodes) 






The second example problem for compressible flow computa- 
tions using 3DBR COMOC corresponds to analysis of turbulent, 
binary three-dimensional mixing in a supersonic boundary layer 
flow field. The analyses include cold flow studies with com- 
parison to data as well as equilibrium combustion. The impetus 
for these studies is the hydrogen fueled scramjet engine, a 
prominent candidate for propulsion of the next generation of 
hypersonic cruise vehicles [Ref. 21, 22]. Many alternative 
designs have been proposed, but all exhibit the commonality that 
fuel introduction arrangements consist of rows of circular, 
choked-flow injector orifices mounted flush or normal to the 
combustor wall or in fins spanning the combustor inlet. The 
pattern of fuel injection, hence three-dimensional mixing, can 
exert significant influence on combustor performance. Figure 17 
illustrates an experimental apparatus used by Rogers to experi- 
mentally probe the three-dimensional cold mixing region down- 
stream of transverse hydrogen injection from a single discrete 
orifice [Ref. 23] , and mul ti pi e laterally-disposed orifices 
[Ref. 24], into a Mach 4 airstream on a flat plate. Detailed 
numerical predictions of turbulent, three-dimensional mixing in 


x i 



Figure 17. Three-Dimensional Flow Field Downstream of 
Transverse Injection from Discrete Orifices 
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this flow field have been performed, using 3DBR COMOC, in a so- 
lution domain spanning up to 120 injector diameters downstream, 
for both the single and multiple jet geometries. The scaled 
discretization of the symmetric half-plane into 100 finite ele- 
ments, Fig. 18, was formed by the automatic discretizer. The 
turbulent boundary layer thickness at the injector station 
(without injection) was equal to 2.7 injector diameters, and 
the Reynolds number was Re x = .6(8)/m. To the first order of 
approximation, these data correspond to isoenergetic binary 
diffusion in a constant pressure flow field. Hence, numerical 
integration was required for u-| and a single species mass frac- 
tion plus the continuity equation for U£. 

Closure of the equation system requires a turbulence model 
for three-dimensional boundary region flow of this type. A 
prototype eddy viscosity model was developed to reflect mass 
flux differences between the main flow and the jet and the tur- 
bulence due to the presence of the wall [Ref. 25]. Directly 
above the mixing region, turbulent dispersion was assumed to 
primarily reflect differences in mass flux. Outside the mixing 
region but near the wall, the turbulence was assumed due solely 
to boundary layer phenomena. Within the mixing region both 
mechanisms are assumed active. Therefore, near the wall where 



Figure 18. Finite Element Discretion of Symmetric Half-Space 
of Single-Jet Injection Geometry 
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mixing length theory is assumed valid, the eddy d i f f u s i v i ty of . 
mass is of the form 

Sc ^ 1 (65) 

u, is the mean longitudinal velocity, Z is the mixing length, 
id is van Driest's damping factor [Ref. 26], and y is the inter- 
mittency factor, empirically modeled in a number of ways [Ref. 
27]. For this study, it was evaluated as 


2 

Z my 


8U- 


3x, 


1 

Y = 

1 + 0 . 01 ? 


( 66 ) 


where r, = X 2 /S 2 » and 63 is the value of X£ where the hydrogen 
mass fraction is one-half its wall value. In the outer region, 
the eddy diffusivity of mass is assumed proportional to the 
mass defect of the form 


c 2 = K Yl/pL (67) 


where K is an empirical constant, y is the i ntermi ttency factor, 
the characteristic length L is defined as. the half height of the 
mixing region on the centerplane (x '3 = 0 ), and I is the three- 
dimensional mass defect evaluated as 


I (*-| ) 


|p U 1 - Poo U J dT 

R 


( 68 ) 


A subroutine was coded for addition to 3DBR-C0M0C, to evaluate 
Eq . (65)-(68) along columns of nodes parallel to the X 2 axis, 
see Fig. 18. The effective viscosity of mass mixing was laminar 
plus selectively e] or e 2 > and the transition from e-j to was 
internally signaled. So as to provide the potential user with 
an example of construction of a fairly complex subroutine 
addition, this code forms a part of the data deck for the second 
standard test case for 3DBR-C0M0C. 

An extensive computational program was conducted using 
3DBR-C0M0C to validate use of the three-dimensional boundary 
region equations for this problem geometry, and to evaluate 
the governing influences on turbulent mixing in the three- 
dimensional region within the constraints of the prototype eddy 
viscosity model. The existence of extensive cold flow data 
[Ref. 23 , 24] for longitudinal velocity, u-i, and hydrogen mass 
fraction, Y H , provides the means to evaluate the accuracy and 
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consistency of the numerical predictions; the complete dis- 
cussions of results is in Ref. 25. Briefly, for the first 
series, the experimental data at 30 diameters downstream of 
injection (i.e.,xi/D = 30) for a dynamic pressure ratio (q r ) 
of unity, were interpolated at the nodes of the finite element 
discretization, Fig. 18, using cubic spline procedures. The 
numerical solution was carried to x-|/D = 60, and the empirical 
constant in Eq. (67) set to K = 0.1 by "best" agreement with ' 
data. As shown in Fig. 19, agreement with the superimposed data 
along the symmetry plane, X 3 = 0, is excellent. Transition from 
mixing length to mass defect occurred between 0.6 and 1.0 injec- 
tor diameters above the plate across the entire pattern. The 
predicted lateral spreading of the hydrogen (parallel to the X 3 
axis) is in good agreement with the data spread near the wall, 
but is underpredicted in the mid-region of the pattern. Transi- 
tion from the initial distribution, and detail on solution ac- 
curacy, are presented in Fig. 20, which compares data to computed 
concentrati on profiles along planes x^ = constant at x 1 / D = 60. 
The underpredicted lateral diffusion is prominent, stemming in 
part from the omission of three-dimensional influence in the 
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Symbol Description 

Initial Conditions, x^/D = 30 

A Flaw Data Measurements 

^ Best Symmetry Plane Determination 

■ ■ COMOC Computation 



Figure 20. Computed Single-Jet Hydrogen Mass Fraction 
Distribution at x-j/D = 60, q r =1.0 

turbulence model. The experimental data show that in this re- 
gion, u i f 2 and Ml , 3 are of equal magnitude, as observed in Fig. 
21 , which is a three-dimensional surface plot of the longitudinal 
velocity distribution at x q / D = 30, as observed from a vantage 
point beneath the plate surface. The superimposed grid coincides 
with the finite element discretization (omitting diagonals), and 
the hydrogen jet is imbedded within the centroidal indentation. 
Obviously, three-dimensional effects are important, and should 
form an integral part of future three-dimensional turbulence 
modeling. Figures 22 and 23 compare the COMOC computed solutions 
for hydrogen mass fraction and longitudinal velocity distribu- 
tion to data at x-|/D = 120, as well as the initial distributions. 
Agreement is generally quite good throughout. 

This computational study was repeated using for initial 
conditions the experimental data for a row of orifices aligned 
perpendicular to the main flow vector with a uniform separation 
distance of 12.5 orifice diameters [Ref. 24]. For this study, 
the finite element discretization of Fig. 18 was simply scaled 
to span the symmetric half zone between ejectors, using the 
automatic discretizer. The vanishing normal gradient boundary 
condition was then applied to both lateral sides of the compu- 
tational region. Shown in Fig. 24 and 25 are the COMOC com- 
puted hydrogen mass fraction profile distributions at stations 
X l/D = 60 and x q / D = 120 compared to data. Figure 26 displays 
the more familiar contour plot at x q / D = 120. These results 
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Transverse Displacement - x^/D 


Centerplane 


Figure 21. Isometric View of Longitudinal Velocity Surface for 
Single-Jet Configuration, x-j/D = 30 


Symbol Description 


Initial Condition, Xj/D = 30 

O Raw Data Measurements 


Best Symmetry Plane Determination 
COMOC Computation 


Xq/'D = 0 


x,/ D = 1.0 


Xo/D = 2.5 


Hydrogen Mass Fraction Y (%) 

Figure 22. Computed Single-Jet Hydrogen Mass Fraction 
Distribution at x-,/0 = 120, q = 1.0 




— — Initial Condition, Spline Fit, x^/D = 30 

Symbol Description g Experimental Data, x-j/D = 120 

— COMOC Computation, x^/D = 120 



Figure 23. Computed. Si ngl e- Jet Longitudinal Velocity 
Distribution at x-j/D = 120, q r = 1.0 

— Symbol 1 Description - - - - - — — 

Initial Conditions, x^/D = 30 ; 

O Raw Data Measurements 
• Best Symmetry Plane Determination 
COMOC Computation 



Figure 24. Computed Multijet Hydrogen Mass Fraction 
Distribution at x-|/D = 60, q ^ = 1.0 
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Initial Conditions, Xj/D = 30 
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Hydrogen Mass Fraction - Y (%) 

Figure 25. Computed Multijet Hydrogen Mass Fraction 
Distribution at x-j/D = 120, = 1.0 
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were obtained using the identical mixing model of the single 
jet, i.e., K = 0.1, Eq . (67), with transition from mixing length 
to mass defect occurring in the region 0.6 £ X 2 /D £ 1.1. Figures 
24 and 25 indicate that centerplane diffusion is somewhat over- 
predicted, while a considerable improvement between computations 
and data has occurred in the lateral region. Figure 26 illus- 
trates how the computed contour patterns merge between jets for 
the multiple injector configuration. Agreement between computed 
and measured velocity distributions at both downstream stations 
was exce l lent. 

Detailed volumes of experimental data, of the tyoe utilized 
to start these discussed solutions, are typically not available 
for complex reacting flow fields. Assuming that the foregoing 
studies have indeed verified the appropriateness of the differ- 
ential equation system, as well as a limited validity for the 
turbulence model, methods for starting a three-dimensional solu- 
tion with combustion are sought. One technique that shows some 
promise, and also admits numerical evaluation of its appropriate- 
ness, is the "virtual source." The theoretical hypothesis is 
simply that the complex transverse injection process can be com- 
putationally replaced by a hydrogen jet imbedded within a boundary 
region flow, and that the distinguishing features of this virtual 
source are solely dependent upon freestream and injectant param- 
eters. The derived model [Ref. 25] captures the essence of the 
barrel shock-Mach disk hypothesis for i njectant-f reestream equi- 
libration [Ref. 28]. To establish computational verification of 
this concept, the cold flow studies were repeated for the virtual 
source established in the plane of injection, i.e., x,/D = 0.0. 

It was assumed to be of elliptical cross-sectional shcipe with the 
minor axis parallel to X 2 » and composed of 100% hydrogen with 
longitudinal momentum determined from the dynamic pressure ratio, 
q r . Computational evaluations of the concept were made for the 
three values of q r for which data exist. Shown in Fig. 27 are 
typical results, obtained for virtual source simulation of the 
multiple injection configuration for q r = 1.0. Superimposed are 
appropriate experimental data [Ref. 23, 24] for the key comparison 
bases of, 1) longitudinal trajectory of maximum hydrogen concen- 
tration, 2) height above the plate of the maximum concentration 
trajectory, and 3) the lateral spreading of the diffusion pattern. 
The peak hydrogen mass fraction is observed to drop precipitously 
from its initial 100% concentration, but to promptly level off in 
good agreement with the multijet data. The deoendent variable 
gradients associated with this solution were quite large, yet 
3DBR-C0M0C maintained stable solution behavior using the automatic 
step-size constraint. The trajectory of the maximum hydrogen con- 
centration above the plate surface is similarly in good agreement 
with multi-jet data. The computed lateral spreading agrees well 
with multi-jet data to x g / D = 30, but is progr es s i vel y underpre- 
dicted (maximum 15%) as the solution continues downstream. Simi- 
lar agreement trends with data were recorded for the "softer" 
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injection case,- q r = 0.5. However, fairly large disparities (up 
to 40%) occurred for the "hard" jet, corresponding to q r = 1.5,. 
but for which the entire .theoretical concept becomes rather sus- 
pect. Figure 27 is completed with a plot of n> an integral 
"mixing efficiency" parameter, [Ref. 24]. This parameter is 
defined as the percentage of molecular hydrogen that would be 
lost to the computed frozen flow H? mass fraction distribution 
if all available H 2 (or 02, depending on the limiting reagent) 
was converted to H 2 O within the computed velocity distribution. 

The virtual source turbulent transport model was identical 
to the prior cold flow studies with two minor exceptions. The 
mixing length hypothesis, e] , was uniformly enforced until the 
minimum velocity in the virtual source depression accelerated to 
within ~2% of the corresponding boundary layer velocity without 
injection. This occurred within 8 diameters downstream of the 
injection Plane for all q r . Downstream of X]/D = 8, transition 
from e] to £2 occurred within one diameter above the Plate sur- 
face. Due to the rather small initial density within the virtual 
source, the initial computed mass defect was d i s oropor t i ona 1 ly 
small. From experimentation, a smaller constant (K = 0.05) was 
found uniformly effective for the three studies. 

With the change of one input flag, the virtual source 
simulation was repeated with equilibrium combustion allowed. 

Shown in Fig. 28 are the typical results of this computational 
simulation, on the multiple comparison bases previously dis- 
cussed for the cold flow tests. The precipitous drop in the 
peak hydrogen mass fraction concentration parallels the cold 
flow experience, but the levels downstream of x-|/D = 7 lie well 
below the non-reacting experience. The trajectory of the peak 
hydrogen level is essentially parallel to the plate surface 
through x-j/D = 30. The lateral spreading of the jet is con- 
siderably more pronounced for X]/D>10 in comparison to the cold 
flow tests. The mixing efficiency parameter, n, reaches 100% 
at x 1 / D = 20, well ahead of the cold flow simulation. These 
differences reflect, in large part, the considerably different 
temperature and density distributions. As shown in Fig. 28, 
ignition occurs immediately downstream of simulated injection, 
and the peak temperature rapidly climbs to 2000°K (~3800°R). 

It remains at this level until stoichiometric mixing is com- 
Dleted, i.e., n = 100%; thereafter, it continuously decreases 
as local heat addition from combustion is unable to balance 
diffusion effects. A typical CPU time on the IBM 360/65 for 
execution of this test with combustion was about 1275 seconds. 
Execution of a cold flow counterpart (to x-j/D = 30) required 
about 950 seconds. The 325 seconds additional for combustion 
reflects the CPU required to execute the temperature iteration 
loop in the combustion subroutine for this particular test case. 
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DATA DECK PREPARATION 


The data decks that generate the nominal Mach 5 isoenergetic 
boundary layer flow, and the three-dimensional virtual source 
simulation, come as standard test cases for 3DBR COMOC. The 
listings of these data decks are included in Appendices A and B. 
Another problem specification can be readily adapted from these 
decks, since approximately one-third of a data deck is associated 
with standard call sequences as well as output format specifica- 
tion and arrangement instructions. These standard data should 
not be altered without reference to the programmer's manual for 
3DBR COMOC [Ref. 1]. The following discussion covers general 
details, and illustrative examples are pertinent to the data deck 
for the Mach 5 test case. Comments and descriptions should be 
interpreted with reference to Appendix A. Subsequently, the 
alterations required to establish the non- un if orm ly discretized 
virtual source problem data deck from the Mach 5 test case are 
presented and discussed. 

Input preparation is subdivided into four phases. 

Phase I . Reference Conditions and Control 
Parameter Specification 


Call 

Parameter 

Code 

Function 

FEBL 



Starts execution of COMOC 

COMTITLE 



Reads one title card to appear on cover 




page of output 

FENAME 



Initialization 

& NAME 01 

NEQKNN 


Integer parameter input 

Number of dependent variables to be in- 




tegrated in XI direction 


IGAS 

0 

Isoenergetic flow with constant Cp 



1 

General flows 


IFR 

0 

Equilibrium composition (IGASeI) 



1 

Frozen composition 


KDUMP 

0 

Suppress debug output 



1 

Print debug output 


NPVSX 


No. of entries in pressure table 


NSCX 

0 

Uniform X3 interval in discretization 



1 

Non-uniform X3 interval in di scret i za t ion 


NSC Y 

0 

Uniform X2 interval in discretization 



1 

Non-uniform X2 interval in di scret i zat i on 

&NAME02 

UINF 


Floating point parameter input 
Reference (freestream) velocity (F/S) 


T0FINF 


Reference stagnation temperature (°R) 


REFL 


Reference length ( F ) 


T0 


Initial XI station (F) 


59 




Function 


Call Parameter Code 


FEDIMN 


LINK! 


TD 

DELP 

EPS 

V START 

XSCALE 
YSCALE j 
CPA.CPH] 
T0A,T0H \ 
XMA.XMH 


Length of XI solution, starting at T0 
(F) 

Percent of TD at which output is desired 

Integration control parameter (.01 to 
.0001 ) 

Percent of TD at which transverse veloc- 
ity ( U 2 ) computation starts 

Multipliers to convert discretization 
to feet 

Specific heats, stagnation temperatures, 
and molecular weights for two-compo- 
nent, i soenerget ic , frozen flow mixing 
(IGAS=0) 

Generate vector lengths and array entry 
points. 


Phase II . Finite Element Discretization 

1 This call generates the finite element discretiza- 
tion of the X2X3 plane. The data are read in free 
format fields. A "T" terminates any sequence. 


A. Automatic Uniform Discretization 
Occurs for NSCX = NSCY = 0 

Set XSCALE = desired element width in the X3 direction 
Set YSCALE = desired element height in the X2 direction 
Read selection keys 
e.g. YSCALE = 0.004 
XSCALE = 0.002 
1 21 , 1 2 , 

T 

Generates discretization made up of 21 node rows x 2 
node columns, or 40 elements (x 1 element). Elements 
are 0.004 F high by 0.002 F wide. 


B. Automatic Non-Uniform Discretization 
Occurs for NSCX = 1 and NSCY =1 

Set X3 discretization first, X2 discretization second. 
Data are used in sets of 3 integers at a time. First 
integer identifies finite element interval concerned, 
next two indicate element width (or height) as ratio 
in feet, e.g., 3 1200 = 3/1 200 (F ) . 
e.g. 1 3 1200,2 1 600,3 5 1200, . . . 

T 

1 1 600,7 1 600,8 7 1200, . . . 

T 

1 11, 1 4, 

T 
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This generates a finite element discretization of 11 
node rows x 4 node columns. The element widths 
(intervals between node columns) are respectively 
3/1200 (F) , 1/600 (F) , .... The height of the first 
7 element rows is uniformly 1/600 (F), eighth is 
7/1200 ( F ) , etc. 

Phase III . Output Specification 

Following the discretization phase, the user can input up to 10 
title cards to head each generated output sequence. 

The next -65 input cards specify output format, see Appendix A, 

and are typically not to be changed without reference to the 
programmer's manual. 

Up to 10 title cards can follow the standard output specification 
to fully describe the problem being solved. This output 
will occur once, directly after printing of the cover page. 

DONE Calls end to output specification phase 

Phase IV . Solution Parameters, Boundary Conditions, 
and Initial Distributions 

Call Function 

VX3ST Establishes NPVSX entries into static 

pressure table as function of XI 
e.g. 11*10.05 0.1 Eleven pressure values at intervals AX1 

of 0.05, starting at XI = 0.1. 

VPVSX Read pressures in PSFA 

e.g. 4.3494 3.41... 

I P I NT -1 Standard Input consisting of integer 

array of numbers corresponding to depen- 
dent variables. Program will integrate 
first NEQKNN of them, plus U2. 

KBNO (N) KBNO ( N ) establishes fixed boundary 

conditions for dependent variables N 
through NN. 

KBNO (NN) 

e.g. KBNO 4 Fixes variable 4 nodes on bottom of 

BOTTOM DONE discretization at their initial input 

values. 
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Call 


Function 


ICALL 

ICALLS 

LINK3 

LINK1 

-1 

-1 ! 
4 ; 

3 j 

! 

i 

> 

i 

Fixed calling sequence for internal 
evaluations, not to be changed. 

VTEMP 

-58 


Read initial nodal total temoerature 
distribution. Non-dimensional ize 
entries by number in location 58 (TREF) 

e.g. 

VTEMP 

82*1800. 

T 

Read: first 82 nodes at T = 1800°R 

0 

VYY 

-(x) 


Reads initial conditions for dependent 
variable N. Non-dimensional i ze entries 
by number in | X | . 


VYYEND (N) 


e.g. 

VYY 

-27 

Initial U2 di str i bution i s all zeros 


42*0.0 

T 

VYYEND 

2 

Non-d imens i onal i ze entries by number 
in location 27 (UREF) 

e.g. 

VYY 

-27 

Initial 111 distribution is zero at 


2*0.0 

2*1 654. . . 

first two nodes, 1654 F/S at second 


T 

72*4004.8 

two,..., last 72 nodes have 4004.8 
F/S. Non-dimens i ona 1 i ze entries by 


VYYEND 

1 

number in location 27 (UREF). 


QKNINT 

DESCRIPT 

DONE 

DESCRIPT 3 


Standard completion of data deck 


DONE 

COMOC 

END 


' EXIT 


If a second test case is desired, 
insert data deck starting with COMTITLE 
before EXIT card. 
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Listed in Table 9 are the changes to the Mach 5 test case 
data deck required to establish the three-dimensional virtual 
source data deck. The complete listing of the latter is included 
as Appendix B. The following explains the alterations with 
respect to input phase and the line numbers in Table 9. 


Input 

P h a s e Line No. Description 


I 


II 

III 



Title card for output cover page 
Reference condition and control parameters 
for combustion calculations using five 
dependent variables 

Form non-uniform discretization, using 11 
node rows * 6 node columns, producing 100 
finite elements 

Title card to head each output call ■ 
Detailed problem description 


IV 


22 ' 

23 Entry locations of longitudinal pressure 

24 distribution (constant ) 

25 66 nodes have uniform stagnation temperature 

26} Initial U1 distribution 


36 / 

37 Initial U2 distribution is zero 

38 Initial 02 distribution reflects location 

ofvirtualsource 

39 Initial N2 distribution 

40 Initial H2 distribution 


Hence* establishing the data deck for a multiple dependent vari- 
able, three-dimensional problem using a non-uniform finite ele- 
ment discretization is readily accomplished. In this case, only 
forty data card changes were required, using the two-dimensional 
Mach 5 data deck as a master deck. 
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Table 9 

DATA DECK CHANGES TO PRODUCE VIRTUAL SOURCE SIMULATION 
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CONCLUDING REMARKS 


This report documents the theoretical foundation and the 
mechanical structure of the Three-Dimensional Boundary Region 
Variant of the COMOC computer program system. A unified effort 
has been made to generate a computational capability that can be 
addressed to a wide range of problems involving comnlex three- 
dimensional flow fields without requiring undue mathematical 
prowess on the part of the user. The success of this type of 
venture can only be measured by the degree to which these goals 
are approached or attained. As with any large computer program, 
it has been debugged to the extent of the specific problems 
already explored. Hopefully, if bugs remain, they are not of 
such a debilitating nature as to severely limit the usefulness 
of the program. In this regard, it is suggested that the poten- 
tial user first experiment with a few problems whose solution 
character is known, in order to attest to the program's perfor- 
mance for a particular problem class. 
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APPENOIX A 


DATA DECK LISTING FOR MACH 5 
TWO-DIMENSIONAL FLOW CHECK CASE 


FEBL 

CGMT ITLE 


CHECK CASE, TWO DIMENSIONAL 

SUPERSONIC FLOW 

WITH PRESSURE GRADIENT 


EENAME 





CNAMEO 1 

NE0KNN=2, 

I GA S=0, 

I FR=0 , 

K DUMP = 0, 


NP VSX= 11, 

N SC X= 1 , 

NSC Y= 1 , 


£ END 
CNAMEO 2 

' U IN F = 40C4. 8 , 

TOF I NF =1 800. , 

REF L= .0 132 



TQ = 0. 1, 

TD=0. 2, 

DEL P=5 .0, 

EPS = . 0 1, 


VSTART=5.0, 

C PA =0. 24 , 

C PH=3 .445, 

T 0A= 1800. 


TQH=C . 0 , 

XMA=2 8. 57 , 

XMH=2 .016, 



X SC ALE = 1 • 0 , 

YSC ALE= 1.0, 



CENO 
FED IMN 
L INK l 

l 



SETUP 

1 2 
T 

1000, 2 2 1CC0, 




1 

1 2 
T 

1000, 21 2 1CCC, 




1 

1 21, 
T 

1 2, 





CHECK CASE, TWO 0 IMFN S I ONA L SUPERSONIC FLOW WITH PRESSURE GRADIENT 
REFERENCE ENC-LISH-FT ENGL I SH-I N M-K-S C-G-S 


CONE 

M PAH A -1 

2 ?. 162 

2 2 2 

2 2 2 

2 2 2 

2 -175 2 

2 2 2 

2 2 2 

2 2 169 

2 2 2 

2 2 2 2 

HOL 1ST 

LENGTH 

VELOC ITY. 

DENSITY 

TEMPERATURE.. .. 

ENTHAL PY 

FROZ.SPEC. HEAT 

VISCOSITY 

LOCAL PRESSURE 


164 
164 
170 
16 5 
2 

176 

177 

168 

2 


163 

163 

174 

2 

2 

2 

178 

167 

2 


FT 

• IN.. 




.M 

. CM . . . 
.CM/S. 
.G/CC. 

FT/S.. .... - 




.M/S 

LBM/FT3. .. . 

. N. A. 



• • • • 

. KG/M3 

RANKINE. .. . 

• N • 



• • • • 

.KELVIN 

. N . A . . 

RTU/LBM. . . . 

. N. A. 



• • • • 

. KJ/KG 

• N • A • • 

btl/lbm-r. . 

. N. A. 



• • • • 

.KJ/KG-K.... 

. N . A . . 

LBM/FT-S. . . 

. N. A. 



• • • • 

. NT-S/M2. ... 

•POISE 

P S F 

. PSI. 



• • • • 

. NT/ M2 

.TORR. 
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LOCAL SOLUTION MACH NC. DPDXl (LBF/FT3 ) MAX. H2 CONC. MIX EFF.IETAI 


Xl/LRFF 


DXl/lREF 

EPSILON 

0X1 M I N/LREF 

ICNUMB 

-l 





2C0 

4*43 

2CC 27 200 2*2 7 



200 

10 

200 

2*10 200 5 8 

2 00 58 

200 

200 

97 

200 

57 2C0 200 30 

200 30 

200 

200 

38 

200 

2*38 



599 






200 

36 

36 

36 36 



200 






61 


100 

134 122 



l 1 

12 

14 

85 



ICSAVE 

- 1 





1248 

285 

320 

284 1 C24 8 



2248 

278 

4248 

5246 8248 



1247 

334 

252 

314 



T U.T, 

FS.RH0.N2, 

V,CP »H TOT «H2 ,0 2 *DI FU 

• PR NO. , LAM. 

VISC..SCT.NO. 

IOMULT 

-I 






14*2 

T U,T,HS,RHC,N2,V,CP t HT0T,H2,02,0IFU,PR NO. »LAM.VI SC. ,SCT .NO. 

RESCRIPT 2 

IJl/UREF T/TREF HSTAT/HREF RHO/PHORFF ELEM.N2 MAS.FRAC 

U 2/UP EF CPF/CPFREF HTO I/HRFF FLEM.H2 M AS . FR ACELEM . 02 MAS.FRAC 

EFF.MU/MUREF EFF,. PRANOTL NO.MU/MUREF EFF. SCHMIDT NO. 

CCMOC 

RESCRIPT 

COMOC CHECK CASE FOR TWO-DIMENSIONAL FLOW WITH PRESSURE GRADIENT. 

A COMPARABLE SIMILARITY SOLUTICN HAS BFEN REPORTED BY CHRISTIAN ET AL, 

ARL 70t,0023. 'SPECIFIC CASE CONSIDERED CORRESPONDS TO MACH NO. 5 BETA=0.5, 

SI 0 )=0 ^'ADIABATIC WALL I. SCLUTIGN STARTED AT X=0.10 FT. WITH SIMILAR 
SOLUTION" iP R'O F I L E . L A !$I'N A’R'SRLQ W WITH VISCOSITY FROM SUTHERLANDS LAW. 

C I SCR Eti/AT ION.' SPANS -THREE TIMES INITIAL BOUNDARY LAYER THICKNESS. 

ISOENERjliVT TC fC.OW-WITH TOTAL TEMPERATURE = 1800 R. 

DONE ‘SiR 

V X 3S;T . 

* ; i, 11*10.05 C.I T XI TABLE FOR PRFSSURE 

■'VPV'SX' 


4.3454 

Y 

3.4 1 

2. 

846 2.46 

2. 2176 

2.02 1.857 

1.73 1.6178 1.53 1.4451 

IPINT 

l 






l 

4 

8 

10 9 

3 2 

. T INTEGRATE 

Ul ,ENTH. ,C2,N2,H2,U3,U2 

KBNO 

BOTTOM 

1 


DONE 





FIXES 

UL 

I VARIABLE 

NO. I) 

ALONG WALL TO 

INITIAL VALUE 

KBNO 

BOTTOM 

2 


DCNF 





FIXES 

U2 

I VARI ARLE 

NO. 2) 

ALONG WALL TO 

INIT IAL VALUE 

KBNO 

BOTTOM 

4 


OONF 
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(VARIABLE NO. 4) ALONG WALL TO INITIAL VALUE 
2 I 1 2 T 

3 T 


12 5 6 


FIXES H 
I CALL -l 

2 5 2 

l CALLS -1 

10 6 4 

LINK3 4 

L INK L 3 

VTEMP -58 

42*1800. 

T INITIAL TOTAL TEMPERATURE PROFILE 
VYY -27 

2*0.0 2*865. 2* 1654. 2*2373. 2*3004. 

2*3982. 2*4CC4. 2 24*4004.8 

T INITIAL Ul PROFILE 
VYY ENO 1 

VYY -27 

2*0.0 2*2.12 2*20.14 2*53.52 2*83.2 

2*253. 2*447. 24*45C. 

T INITIAL U2 PROFILE 
VYY END 2 

VYY 


2*3550. 2*3879. 


2*109.3 2*165. 


42*. 2330 

T INITIAL 02 MASS FRACTION PROFILE 
VYY END 8 

VYY 


42*. 7670 

T INITIAL N 2 MASS FRACTION PROFILE 
VYY ENO 10 

VYY 


42*0. C 

T INITIAL H2 MASS FRACTION PROFILE 
VYY END 9 

QKNINT 
OESCR IP T 
OCNE 


CESCRIPT 3 

REFERENCE LENGTH, LREF 
REFERENCE VI SCOSITY. LAMINAR VALUE 
FVALUATEO AT RFF . TEMPERATURE. 
FREESTREAM VELOCITY AT XO(=UREF! 
STAGNATION TEMPERATURE (CON S TA NT ,=TREF » 
FREESTRFAM DENSITY AT XO(=RHOREF) 
FRFESTRFAM MACH NUMBER AT XO 
STATIC PRESSURE AT XC 
NUMBER OF NODES 
NUMBER OF FINITE ELEMENTS 
DONE 


CGMOC 
ENO 
EX IT 


43 FT. 

38 LBM/FT-S 
27 FT/S 
58 DEG R 
10 LBM/FT3 
154 
9 PS F 
-16 
-14 


DIMEN 

GEOMFL 
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APPENDIX. 6 


DATA DECK LISTING FOR VIRTUAL SOURCE 
THREE-DIMENSIONAL CHECK CASE 


PERL 

COMT ITL E 

COMOC CHECK CASF FOP THREE DIMENSIONAL REACTING BOUNDARY REGION FLOW 


FENAME 

CNAMEOi 

NEQKNN=5, 

IGAS=1 , 

I FR = 0 , 

K DUMP=0» 


NPVSX=2 , 

N SC X= 1 , 

NSC Y= l , 


NEIE2 

CENC 

SNAME02 

= 1, 

U INF =22 72. , 

T CF I NF= 533. 0 , 

REFL=. 003333333, 



TC = 0 .0, 

TD = 0. 10, 

DEL P= 5.0, 

EPS=0.0l, 

CENO 
FEDTMN 
L INK l 

XSCALE=C. 003333333, 

1 

Y SC ALE=0. 003333333 

, VST ART= 10 1 .0 , 

SETUP 

1 75 

100, 2 50 100. 3 

125 IOC, A 150 

100, 5 225 100, 



T INCREMENTS BETWEEN X3, NODE-NUMERATOR-DENOMINATOR 

I 5 IC, 7 5 IC, 8 125 100, 9 175 100, 10 250 100, 

T INCREMENTS BETWEEN X2 

1 11,1 t , 

T 11 ROWS AND 6 COLUMNS NORMALIZED BY LREF, HENCE X-Y SCALES = LREF 

CHECK CASE, THREE DIMENSIONAL REACTING BOUNDARY REGION - VIRTUAL SOURCE 

REFERENCE ENCLISH-FT ENGLISH-IN M-K-S C-G-S 


■DONE 

M PARA - 1 

2 7 

2 2 

2 2 

2 2 

2 -175 

2 2 

2 2 

2 2 

2 2 

2 2 


162 1 6 A 163 

2 16 A 163 

2 170 17A 

2 165 2 

2 N 2 2 

2 176 2 

2 177 178 

169 168 167 

2 2 2 

2 2 


MOL 1ST 


LENGTH 

VELOC I TY 

• FT 

.FT/S. ... 

. IN.. . 


.M..... 

. M/S 

.CM. . . 
.CM/S . 
.G/CC. 

OENSITV 

•LRM/FT3. .. . 

. N. A... 


• KG/M3 

TEMPER ATURE.. .. 

.RANKINE... . 

• K« A* • 


.KELVIN 

.N. A. . 

ENTHALPV 

.8TU/LBM. .. . 

. N. A. . 


. KJ/KG 

. N. A . . 

FROZ.SPEC. HEAT 

.8TL/LBM-R. . 

. N. A. . 


.KJ/KG-K. ... 

.N. A. . 

VISCOSITY. ..... 

• LBM/FT-S.. • 

.N. A. . 


.NT-S/M2. ... 

.POISE 

LOCAL PRESSURE 

.PSE 

. PSI. . 


. NT/ M2 

.TORR. 

LOCAL SOLUTION 

MACH NO. 

OPDXl ILBF/FT3) 

MAX. H2 CONC. 

MIX EFF 


I ETA ) 


71 



Xl/LREF 


OXl/LREF EPSILON DX1HIN/LPEF 

TONUMB 

-l 



200 

4*43 

200 27 200 2*27 

200 

10 

200 

2*10 200 58 200 58 200 

200 

97 

200 

97 200 200 30 200 30 200 

200 

38 

200 

2*38 

999 




200 

36 

36 

3 6 36 

200 




61 


100 

134 122 

1 1 

12 

14 

E 5 

IOSAVE 

- 1 



1248 

285 

320 

284 10248 

2248 

278 

4248 

9246 8248 

1247 

334 

292 

314 

T U.T, 

HS , R HO , N2, 

V,CP ,HT0T,H2,02,0IFU,PR NO. , L A M. V I SC . , SCT . NO . 

IOMULT 

- 1 




14*2 

r U, T, HS,RH0,N2, V,CP,HT0T,H2,02,0IFU,PR NO. , LA M. V I SC . , SCT .NO . 

OESCRIPT 2 

Ul/UREE T/TREF H ST AT /HREF RHO/RHOPEF ELEM.N2 MAS. FRAC 

0 2/ UR EE CPF/f. PFREF HTOT/HREF ELEM.H2 M AS . F R ACELEM . 02 MAS. FRAC 

EFF.MU/MtJREF EFF. PR AND TL NO. MU/ MUREF EFF. SCHMICT NO. 

COMOC 
OESCR IPT 

CHECK CASE, THREE DIMENSIONAL REACTING BOUNDARY REGION - VIRTUAL SOURCE 
(H2/02/AIR SYSTEM WITH EQUILIBRIUM CHEMISTRY). PROBLEM CONSIDERED 
REPRESENTS TRANSVERSE H2 INJECTION INTO A SUPERSONIC AIR STREAM 
CHARACTERISTIC OF SCRAMJET FUEL I NJECTI CN, SEE ROGERS NASA TNC-6114, 

1971 AND NASA TNO-6476, 1971 FOR EXPERIMENTAL STUDY OF THIS PR08LEM. 

TURBULENCE MODEL EMPLOYED IS DESCRIBED IN USER»S MANUAL NASA CR-132450, 1974. 
CALCULATIONS ARE STARTEO USING VIRTUAL SOURCE CONCEPT TC REPLACE 
COMPLEX NEAR INJECTION FLOW FIELO. 

DONE 

VX3ST 


— 

0.0 

100. 

T XI TABLE 

FOR 

PRESSURE 





VPVSX 








. 


193. 

193. 

T PRESSURE 

TABLE PSF 





TP INT 

-1 









l 

4 8 

10 9 

3 

2 

T INTEGRATE 

Ul » E NTH 

. ,02*N2,H2,U3,U2 


K PNO 
BOTTOM 

1 

DONE 









FIXES Ul 

(VARIABLE 

NO. 

1) 

ALONG WALL 

TO 

I N IT I AL 

VALUE 


K BNO 
BOTTOM 

2 

DONE 









FIXES U2 

(VARIABLE 

NO. 

2) 

ALCNG WALL 

TO 

INITIAL 

VALUE 


K BNO 
BOTTOM 

4 

3DCNE 









FIXES H 

(VARIABLE 

NO. 

41 

ALONG WALL 

TO 

I N IT I AL 

VALUE 


72 



KRNO 

BOTTOM 


8 


30CNE 



FIXES H 

(VARIABLE NO. 

8) 

ALONG 

WALL 

TO 

INIT I AL 

VALUE 

KPNO 

9 








ROT TOM 

• FIXES H 

3DCNE 

(VARIABLE NO. 

9) 

AL CNG 

WALL 

TO 

INITIAL 

VALUE 

KBNO 

BOTTOM 

10 

3D0NE 








F I X c S H 

(VARIABLE NO. 

10 

ALONG 

WALL 

TO 

INITIAL 

VALUE 

ICALL 

-l 








2 

5 2 

2 11 2 T 







ICALLS 

- 1 








lb 

6 A 

12 5 6 3 T 








L INK 3 A 

L INK l 3 

VTEMP -58 

66*533. 

T INITIAL TOTAL TEMPERATURE PROFILE 
V YY -27 


6 * 0.0 

6*1503. 

6*1660. 

2*1550. A* 1759. 

2*1550. A* l 83 3 . 

2*1550. A* 1 852 . 

2*2272. A* 1 5A 2 . 

2*2272. A* 1585. 

2*2272. A*2C7A . 

2*2272. A* 2 165. 

6*2272. 

T INITIAL U1 PROFILE 
VYYEND 1 

VYY -27 


66 * 0.0 

T INITIAL U2 PROFILE 
V vy em n 2 

VYY 

19*. 233 2*0.0 A*. 233 2*0.0 A*. 233 2*0.0 

T INITIAL 02 MASS FRACTION PROFILE 
VYYEND 8 
VYY 


18*. 767 
T INITIAL N 2 


2*0.0 A*. 767 

MASS FRACTION PROFILE 


2 * 0.0 


A*. 767 2*0.0 


VYYEND 10 
VYY 


19*0.0 2*1.0 A*0. 0 2*1.0 

T INITIAL H 2 MASS FRACTION PROFILE 
VYYEND 9 
QKN INT 
OESCRIPT 
DONE 

DESCRIPT 3 

REFERENCE LENGTH, LREF 

REFER P NCE V ISCO S ITY . LA M l NAR VALUE 

EVALUATED AT REF. TEMPERATURE. 

FREESTREAM VELOCITY AT XC(=UREF) 

STAGNATION TEMPERATURE ( CON S TA NT ,= TREF ) 
FREESTREAM DENSITY AT XO(=RHCREFI 
FREESTREAM MACH NUMBER AT XO 
STATIC PPESSURE AT XC 
NUMBER OF NODES 
NUMBER OF FINITE ELEMENTS 
DONE 
END 


A*0 .0 2*1.0 


A3 

FT . 

38 

LBM/FT-S 

27 

FT/S 

58 

DEG R 

10 

LBM/ FT 3 

15A 


9 

PS F 

-16 


-1A 



EXIT 


DIMEN 

GEOMFL 


3A* .233 


3A*.767 


3A*0 .0 
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